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CAVITATION IN LIQUID CRYOGENS 
I - VENTURI 

J. Hord, L. M. Anderson, and W. J. Hall 
1. SUMMARY 

This document constitutes the first of four volumes to be issued 
on the results of continuing cavitation studies. It is an extension of a 
previous study that defined the cavitation characteristics of liquid hydro- 
gen and liquid nitrogen flowing in a transparent plastic venturi. Thermo- 
dynamic data, consisting of pressure and temperature measurements 
within fully developed hydrogen cavities, are reported here. In the pre- 
vious study, it was concluded that the measured temperatures and pres- 
sures within the central and aft regions of the cavities were generally 
not in stable thermodynamic equilibrium. Because these results were 
not anticipated, it was decided that additional, more precise tests were 
in order. Accordingly, the plastic venturi has been equipped with addi- 
tional pressure and temperature sensing ports, the temperature instru- 
mentation improved, and additional tests performed. The new data posi- 
tively confirm the older data. The cavity pressure depressions (bulk- 
stream vapor pressure less measured cavity pressure) increased with 
increasing velocities, cavity length, and fluid temperature. Minimum 

measured cavity pressure was less than bulkstre-am vapor pressure by 

/ 2 

as much as 15. 13 psi (10. 44 N/cm ); measured temperatures and pres- 
sures within the vaporous hydrogen cavities substantiated thermodynamic 
metastability in the central regions and trailing edges. 

Existing correlative theory is used to obtain equations that corre- 
late the new (and old) experimental data for developed cavitation in liquid 
hydrogen. The new equations are shown to be compatible with the older 
data and with the work of others. 
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Details of the test apparatus, test model, instrumentation, test 
procedure, data analysis, and correlative techniques are discussed. 
Experimental data resulting from this study are presented in tabular 
form over the range of experimental temperatures, 35. 77 to 41. 36 R 
(19. 87 to 22. 98 K) and inlet velocities, 98. 8 to 194. 0 ft/ s (30, 1 to 
59. 1 m/s). 

2. INTRODUCTION 

Vaporous cavitation is the formation of the vapor phase within a 
flowing liquid, due to a reduction in pressure. Since the formation and 
collapse of vapor cavities alters flow patterns, cavitation may reduce 
the efficiency of pumping machinery [1]^ and reduce the precision of 
flow measuring devices. Collapse of these vapor cavities can also cause 
serious erosion damage [2] to fluid-handling equipment. While the non- 
cavitating performance of hydraulic equipment may be predicted from 
established similarity laws, cavitating performance is much more dif- 
ficult to predict from fluid -to -fluid. Recent advances in this area have 
been made by NASA-JLeRC personnel [3-5], but additional work is re- 
quired to improve the current technique for predicting cavitating perform- 
ance of equipment from fluid -to -fluid. The effects of fluid properties on 
cavitation performance are well recognized [6-15] and require more 
understanding to develop improved similarity relations [15] for equip- 
ment behavior. Considerably more knowledge is needed to extend this 
predictive capability from one piece of equipment to another, i. e. , a 
more general predictive technique, applicable to equipment design, must 
include the effects of equipment geometry and size, in addition to fluid 
properties. 


1 

Numbers in brackets indicate references at the end of this paper. 
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NASA has undertaken a program [1] to determine the cavitation 
characteristics of various hydrodynamic bodies and the thermodynamic 
behavior of different fluids, in an effort to obtain improved design cri- 
teria to aid in the prediction of cavitating pump performance. The experi- 
mental study described herein was conducted in support of this program. 

In the original work [16], liquid hydrogen and liquid nitrogen were 
chosen as test fluids for the following reasons: (1) the ultimate goal of 
this program is to acquire sufficient knowledge to permit intelligent 
design of pumps for near-boiling liquids, and (2) predictive analyses 
indicated [1] that the physical properties of hydrogen and nitrogen make 
them particularly desirable test fluids. The objective of that study [16] 
was to determine the flow and thermodynamic conditions required to 
induce incipient and developed cavitation on the walls of a transparent 
plastic venturi, using liquid hydrogen and liquid nitrogen. The shape of 
the venturi was chosen to duplicate the test section used by NASA [15]. 

One of the most interesting results of the initial work [16] was the 
indication that thermodynamic metastability exists in the midregions of 
the vaporous hydrogen cavities. In those tests, the uncertainty of the 
cavity temperature measurements precluded accurate definition of the 
magnitude of metastability within the cavities. Although the earlier 
work [16] was in itself complete, the instrumentation was less than opti- 
mum. Consequently, it was necessary to increase the number of pres- 
sure and temperature sensing ports to more clearly define the pressure 
and temperature profiles within the cavities; also, recent advances in 
cryogenic thermocouple thermometry and high-speed data acquisition 
systems made it possible to rapidly and precisely measure temperature 
at several locations in the cavity. The plastic venturi has been equipped 
with additional sensing ports, the temperature instrumentation improved, 
and the hydrogen tests repeated; the results of these improved tests are 
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presented herein and they positively confirm the initial data. Pressure 
and temperature profiles within fully developed hydrogen cavities were 
measured and are referred to herein as developed cavitation data. Ven- 
turi inlet velocities were varied from 98. 8 to 194 ft/s (30. 1 to 59. 1 m/s) 
and inlet temperatures ranged from 35. 77 to 41. 36 R (19. 87 to 22. 98 K). 
The bulkstream vapor pressure exceeds the measured cavity pressure 
and the saturation pressure corresponding to the measured cavity tem- 
perature; therefore, the measured pre s sure depressions and the pres- 
sure depressions corresponding to the measured temperature depressions 
within the cavity are called "pressure depressions. " Alternatively, the 
pressure depression may be expressed in terms of its equivalent equili- 
brium "temperature depression. " 

A similarity equation has been developed [15] for correlating cavi- 
tation data for a particular test item from fluid -to -fluid; this correlation 
is also useful in extending the velocity and temperature range of data for 
any given fluid. The experimental data from this study have been used 
to evaluate the exponents on various terms in this correlating equation. 
These new data have also been correlated with the old data [16] using 
the similarity equation. All data reported here are intended to supple- 
ment that given in reference [16] and are tabulated in appendix A. 

Test apparatus, test procedure, instrumentation, and data analysis 
are described in this paper for the sake of completeness and because 
they have been substantially improved since the initial work was reported 
[ 16 ]. 

3. EXPERIMENTAL APPARATUS 
The facility used for this study consisted of a blow-down system 
with the test section located between the supply and receiver dewars; 
see figure 3. 1. Dewars and piping were vacuum shielded to minimize 
heat transfer to the test fluid. Flow control was attained with a motorized 
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Pressure Control 
Valve 








throttle valve and by regulating the supply and receiver dewar pressures 
Pressure and volume capacities of the supply and receiver vessels are 
indicated on figure 3. 1. The receiver dewar pressure control valving 
limited the inlet velocity, V^, to around 200 ft/s (61 m/s) in these hydro 
gen tests. Valves located on each side of the test section permit flow 
stoppage in the event of venturi failure while testing with liquid hydrogen 
A plenum chamber was installed upstream of the test section to assure 
uniform non-cavitating flow at the test section inlet. The supply dewar 
was equipped with a 5 kW heater which was used to heat the test fluid. 

3. 1 Test Section 

A photograph of the test section, as viewed through one of the win- 
dows in the vacuum jacket, is shown in figure 3. 2; this photograph was 
taken during the initial test series [16] and does not show the additional 
pressure and temperature sensing ports. The transparent plastic ven- 
turi was flanged into the apparatus using high compression elastomeric 
n O n rings. Test section details are given in figures 3. 3 and 3. 4. All 
but two of the sensing ports detailed in figure 3. 3 were used during the 
developed-cavity tests to determine the temperature and pressure de- 
pressions within the cavities. Sensor ports 3 A and 3B were not used 
to avoid possible breakage of the plastic. Cavity length was determined 
from scribe marks on the plastic venturi; see figure 3. 2. The theo- 
retical and as-built venturi contours are shown on figure 3.4. The test 
section dimensions were checked by using the plastic venturi as a mold 
for a dental plaster plug. The plug was then removed and measured. 
Pressure distribution for non-cavitating flow across the quarter round 
contour [16 - 18] is shown in figure 3. 5. This pressure profile has been 
confirmed using several test fluids [15], and applies when 

(Re) D MxlO 5 . 
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Figure 3. 2 Photograph of plastic venturi test section installed in system. Note 
counter --used to correlate flow data with film event. 
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Figure 3. 5 Pressure distribution through test section for non-cavitating flow, 










3. 2 Instrumentation 


Location of the essential instrumentation is shown on figures 3. 1 
and 3. 3. All of the error statements given in this section are based on 
an estimated systematic error and an estimated random error. The 
random error, or imprecision, has been assigned a 3 o confidence limit, 
i. e, , the random error cited will include 99. 73 percent of the experi- 
mental observations. 

Liquid level in the supply dewar was sensed with a ten-point car- 
bon resistor rake. Test fluid temperature in the supply dewar was 
determined by two platinum resistance thermometers, see figure 3. 1. 
Fluid temperature at the flowmeter and test section inlet were also 
measured with platinum resistance thermometers. These platinum 
thermometers were calibrated against a secondary thermometer stand- 
ard and were powered with a current source that did not vary more than 
0. 01 percent. Voltage drop across the thermometers was recorded on 
a 4 digit electronic voltmeter data acquisition system. The overall 
uncertainty of the platinum thermometer temperature measurement is 
estimated to be within ± 0. 09 R (+ 0, 05 K), with an allowance of 
± 0. 06 R(± . 03K) for systematic error and ± 0. 03 R(± . 02 K) for 
imprecision. 

Chromel-gold (0. 07 atomic percent iron) thermocouples were 
used to determine the temperature profile within the cavities during the 
tests. These thermocouples more than double the signal voltage of the 
Chromel-constantan wires used in the initial study [16]. The reference 
thermocouples were placed in the plenum chamber beside the platinum 
thermometer used to determine bulkstream temperature at the test 
section inlet. The thermocouples had exposed soldered junctions and 
were constructed from 36 AWG wire to ensure rapid response. The 
detailed construction of these thermocouples is given in appendix B. 


Tradename --See footnote on page 58. 



The thermocouple leads extending from the reference to cavity thermo- 
couples were thermally lagged to the cold pipe and radiation shielded 
with multilayer aluminum foil. The signal leads which penetrated the 
vacuum barrier were also tempered to the cold pipe and radiation 
shielded to minimize heat transfer to the low temperature junctions. 
Details of the thermocouple installation are shown on figures 3. 6 and 
3. 7. The thermocouple circuits were calibrated in situ, over the range 
of experimental velocities and temperatures, from tests involving non- 
cavitating flow through the venturi; i. e. , pre and post calibrations were 
obtained during each developed-cavity test by causing non-cavitating 
flow to occur with only a slight variation in flow conditions. Overall 
uncertainty of the cavity temperature measurements is estimated at 
± 0. 36 R (± 0. 20 K), with an allowance of ± 0. 18 R(± 0. 10 K) for sys- 
tematic error and ± 0. 18R (± 0. 10 K) for imprecision. 

System gage and differential pressure measurements were obtained 
with pressure transducers mounted in a temperature stabilized box near 
the test section. Differential pressure measurements were used where 
possible to provide maximum resolution. The pressure transducers 
were calibrated via laboratory test gages and manometers at frequent 
intervals during the test series. These transducers have a repeatability 
of better than ± 0. 25 percent, and their output was recorded on a mag- 
netic tape data acquisition system with better than ± 0. 25 percent resolu- 
tion. The overall uncertainty of the pressure measurements, including 

calibration and read-out errors, is estimated to be within ± 1. 0 psi 

2 2 
(± 0. 69 N/cm ), with an allowance of ± 0. 2 psi (± 0. 14 N/cm ) for sys- 

/ 2 

tematic error and ± 0. 8 psi (±0. 55 N/cm ) for imprecision. Bourdon 
gauges were used to monitor the tests. 

Volumetric and mass flow rates were determined via a Herschel 
venturi flow meter designed to ASME Standards [19]. The internal 


12 



WALL OF 
PL ASTI C 
VENTURI 



couple Junction 
with wall ) 


>uple 


Shields are electrically 
isolated from each 
other at this point, 



locouple 


s of thermocouples used 
res. 

1 3 




7ZZZZZZZZZ2 ZZ Z2. TZZ3. ZZZZZZZZZ ZZZZZZZZZZZ ZZZZZZZ ^ 



14 


r 


contour of this meter was verified in the same manner as the test ven- 
turi. An error analysis of this flow device and associated pressure and 
temperature measurements indicates an overall uncertainty in mass flow 
of ± 1.0 percent, with an allowance of ± 0. 2 percent for systematic error 
and ± 0. 8 percent for imprecision. The precision of the computed mass 
flow was periodically verified by comparison with the rate of efflux of 
liquid from the supply vessel; rate of efflux at low and intermediate flow 
rates was determined from liquid level-time measurements in the liquid 
level-volume- calibrated supply vessel. 

An electronic pulsing circuit provided a common time base for 
correlating data between oscillograph (used for monitoring of tests), 
magnetic tape data acquisition system, and movie film. The electronic 
pulsing circuit was triggered by the scanner on the multi-channel data 
acquisition system. This 16 channel, magnetic tape recorder system 
was equipped with a multiplexor so that 3Z data channels were sampled 
each second. Calibration information and data on the magnetic tape 
were subsequently used in a digital computer program to provide data 
for analysis. The data were reduced by first viewing film of the test 
event. A solenoid- actuated counter and a small light bulb, installed 
adjacent to the test section--see figure 3. Z--were energized by the 
electronic pulser and appear on the film. Thus, the data retrieved 
from the magnetic tape are reduced at the desired instant of time by 
simply matching the number of voltage pulses which have elapsed. 

3. 3 Visual and Photographic Aids 

Use of a plastic test section, liquid hydrogen, and relatively high 
pressures precluded direct visual observation; therefore, closed-circuit 
television was used to observe the tests. 

Movies of cavitation tests were taken at approximately 20 frames 
per second on 1 6 mm film. The variable speed camera was equipped 
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with a 75 mm lens and synchronized with a high intensity stroboscope, 
providing a 3 microsecond exposure. The stroboscope was situated 
directly opposite the camera lens and illuminated the test section through 
a plastic diffuser mask; this technique provided bright field illumination 
or a back-lighting effect and excellent contrast between vapor and liquid 
phases in the test section. A low-intensity flood light was also used to 
provide continuous back -lighting for television reception. 

4. TEST PROCEDURE 

The following procedure refers to figure 3. 1. The supply dewar 
was filled with test liquid and then some of the liquid extracted through 
valves A and B to slowly cool the test section and piping. Approximately 
2 hours were required to cool the plastic test section without breakage. 
Cooldown was monitored with the platinum resistance thermometer in 
the plenum chamber. Upon completion of cooldown, the remaining con- 
tents of the supply dewar were transferred through the test section into 
the receiver dewar, and then back into the supply dewar again. This 
operation cooled the entire flow system in preparation for a test. The 
test section and connecting piping were kept full of liquid (at near- 
ambient pressure) during preparatory and calibration periods between 
tests; therefore, the plastic venturi was generally colder than the test 
liquid in the supply dewar when a test was started. Next, the liquid in 
the supply dewar was heated to the desired temperature. Supply and 
receiver dewars were then pressurized to appropriate levels, throttle 
valve D was positioned, and flow started by opening valve C. In the 
case of non- cavitating flow, inception was induced by further opening 
valve D and thus increasing the flow velocity until vapor appeared. To 
obtain desinent cavitation or non- cavitating flow from fully developed 
cavitating flow, valve D was further closed until the vapor cavity was 
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barely visible. For developed-cavity tests, valve D was adjusted to 
obtain the desired cavity length. Flow was terminated by closing valve 
C. The supply dewar was then vented and the test liquid transferred 
back into the supply dewar for another test. As previously mentioned, 
the entire test event was recorded on movie film which was subsequently 
used to determine incipience, desinence, and desired cavity lengths. 

A concrete protective barrier separated operating personnel from 
the test apparatus, and all tests were remotely controlled. Valve D was 
electrically driven and the receiver dewar pressure was remotely con- 
trolled by means of a pneumatic transmitter, differential controller, and 
vent valve arrangement, figure 3. 1. The supply dewar and valve C were 
also remotely controlled with pneumatics. 

5. DATA ANALYSIS 

The developed cavitation data for liquid hydrogen are given in 
complete detail in appendix A. Typical profiles of measured pre s sure 
depression are given on figures 5. 1 to 5. 5. Photographs of fully devel- 
oped vaporous cavities in liquid hydrogen are shown on figure 5. 6. Inlet 
velocity and temperature were observed to have very little effect on the 
appearance of cavitating hydrogen; i. e. , the cavities were generally well 
defined and uniformly developed. 

5. 1 Metastable Trends in Vaporous Hydrogen Cavities 

In figures 5. 1 to 5. 5, the data points representing cavity pressure 
measurements have been connected with a smooth curve to facilitate 
comparison with the data points obtained from the cavity temperature 
measurements. The pressure depressions obtained from the cavity 
temperature measurements are, for the most part, greater than those 
derived from the measured cavity pressures. The new hydrogen data 
indicate that the cavity pressures and temperatures are nearly in 
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Figure 5. 1 Pressure and temperature depressions within cavity- 
in liquid hydrogen. 
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Figure 5. 2 Pressure and temperature depressions within cavity 
in liquid hydrogen. 
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Figure 5. 3 Pressure and temperature depressions within cavity- 
in liquid hydrogen. 
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Figure 5. 4 Pressure and temperature depressions within cavity- 
in liquid hydrogen. 
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Figure 5. 5 Pressure and temperature depressions within cavity 
in liquid hydrogen. 
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(. 1): Typical incipient cavitation. 

Note se-ribe marks used to 
identify nominal cavity length. 


(.2): Nominal cavity length, 1. 25 -inch 
(3. 2 cm); V Q = 33. 8 m/s, 

T = 20. 3K. P = 16. 1 N/cm 2 , 
K° = 1. 46. ° 



(. 3): Nominal cavity length, 2-inch 

(5. 1 cm); V Q = 47. 2 m/s, 

T = 20. 4 K, P Q = 25. 1 N/cm, 
K° = 1. 83. 



(.4): Nominal cavity length, 3.25-inch 

(8. 25 cm); V 0 = 62. 5 m/s, 

T = 22. 7 K, P = 41. 0 N/cm 

o # O 9 

K v = 1. 60. 


Figure 5. 6 Photographs showing typical appearance of 
developed cavities in liquid hydrogen. 


23 


II 





AXIAL DISTANCE FROM MINIMUM PRESSURE LOCATION , X , cm 






thermodynamic equilibrium near the leading edge of the cavity where 
vaporization occurs; these data also substantiate metastability in the 
central and aft regions of the vaporous hydrogen cavities. This is an 
important observation, because it was shown in reference [20] that the 
B-factor theory requires only that thermodynamic equilibrium prevail 
during the vaporization proces s; consequently, only the leading edge 
of the developed cavity must be in thermodynamic equilibrium. The 
central and trailing regions of the cavity, where condensation occurs, 
apparently are not in stable thermodynamic equilibrium. 

The magnitude of thermodynamic metastability within the hydro- 
gen cavities is plotted, for all of the new data, on figure 5. 7. The maxi- 
mum overall uncertainty in the data is clearly indicated on figure 5. 7; 
this maximum uncertainty is obtained by the simple addition of the over- 
all uncertainties in P and P . From section 3. 2 of this paper, the 

n n , T ^ 

estimated overall uncertainty in P^ is ± 1. 0 psi (± 0. 69 N/cm ) and the 

uncertainty in T is ± 0. 36 R (± 0. 20 K). Converting the uncertainty in 
n 

T into uncertainty in P^ ^ merely requires evaluation of the slope of 

the saturation liquid-vapor pressure-temperature equilibrium curve 

at the appropriate T^ t The product of this slope and the uncertainty in 

T is the estimated uncertainty in P . Because the uncertainty in 
n n, T 

P varies with temperature level, the data have been separated on 
n , T 

figure 5. 7 into two discrete groups; each group identifies with a pre- 
scribed fluid temperature bandwidth and the associated maximum uncer- 
tainty in the data. The estimated overall uncertainty in P at the 

n , T 2 

lower temperature level is (2. 39 psi/R)x 0. 36 R = 0. 86 psi (0. 59 N/cm ), 

and at the higher temperature level is (3. 7 5 psi/R)x 0. 36R = 1. 35 psi 
/ 2 

(0. 9 3 N/cm ). Errors in the hydrogen property data will be taken as 

negligible in this analysis. The uncertainty bands indicated on figures 

/ , 2 

5. 7 are obtained by adding 1. 0 psi (0. 69 N/cm ) to the uncertainties 

in P listed above. These limits appear quite realistic as they encompass 
n, T 
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practically all of the data at instrument stations 2 and 4, where equili- 
brium apparently exists; the distribution of data points at these two 
stations lends credence to this statement. Stations 6 to 9 exhibit a 
definite trend toward metastable vapor, with the maximum consistent 
metastability occurring near station 6. 

These results confirm the initial work [16] where cavity tempera- 
tures were measured at stations 2, 6, and 8 only. It is felt that the data 
plotted on figure 5. 7 is positive evidence that metastable vapor exists 
in the central and aft regions of the hydrogen cavities. In this plot, no 
attempt was made to separate the effects of velocity and cavity length; 
considerably more data would probably be required, and such effort 
does not appear warranted at this time. From figure 5. 7, it appears 
that fluid temperature does not have a significant effect on the magni- 
tude of P - P _ . 

n n, T 

From figure 5. 7, it can be observed that, within data accuracy, 
thermodynamic equilibrium exists in the frontal regions of the cavity. 
Over the remaining cavity length, it appears that condensation occurs-- 
an important consideration for future modeling theories. Because of 
the importance of this result, one is tempted to hypothesize an explana- 
tion for the results shown on figure 5. 7. Accordingly, the following 
discussion should be considered as conjectural. 

5. 2 Simplified Analysis of Metastable Phenomena 

Flow patterns, circulation, turbulence, and vapor velocity within 
the cavities are unknown; consequently, it is difficult- -if not impossible-- 
to perform a thermodynamic fluid flow analysis on the liquid-vapor ex- 
change. Thus, we shall limit our discussion to a highly simplified analy- 
sis in order to obtain some crude explanatory results. We will assume 
that the developed cavity can be represented by an agglomeration 
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of single spherical bubbles moving with the bulkstream liquid; we can 
then use the existing solutions [21] for bubble dynamics. These solutions 
are available in simplified form, appropriate for the approximate solu- 
tions we seek. 

The following assumptions were used to obtain the simplified ex- 
pressions for vapor-filled bubbles; The flow is taken as steady and 
irrotational; compressibility of the liquid, heat flow, viscosity and sur- 
face tension effects are neglected. Also, a constant pressure and tem- 
perature in the bulkstream liquid is assumed. Inertial forces are con- 
sidered dominant, and thermal boundary conditions and pressure fluc- 
tuations are ignored; for our problem, this constitutes the only serious 
breach between the simplifying assumptions and our intended application. 
We know that heat transfer and pressure gradients play an important role 
in the formation and sustenance of a developed vapor cavity. A more 
rigorous solution, including pressure variations and heat transfer to a 
single bubb le, is a complex task [22] and for the purpose of this dis- 
cussion is not warranted; particularly, since our problem involves an 
annular developed vapor cavity. In simplified form, the cavity growth 
rate approaches the value given by 



and the collapse rate is approximated by 


(5. 2-1) 


R 


C 




(5.2-2) 
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The pressure differences used to calculate cavity growth and 

collapse rates in the foregoing equations have been chosen to take 

advantage of measurements provided by this study. Also, extension 

of the spherical bubble formulae [21] to our annular cavity problem 

requires some improvisation. The conventional bubble formulae [21] 

use P - P, in eq (5. 2-1) and P, - P in eq (5. 2-2). That we 
v,x -t,x v,x v,x 

cannot blindly use these quantities is apparent because 1) P. is 

X 

unknown- -it is not measured and is difficult to calculate because it 

varies continuously with cavity length and is also influenced by the 

unknown cavity thickness, 2) P is not constant--it varies continu- 

v, x 

ously along the annular cavity interface as a result of vaporization (or 

condensation) and changing bulkstream liquid pressure. Extrapolation 

of the bubble formulae to the annular cavity may be accomplished by 

assuming 1) that the liquid pressure, P , is adequately approximated 

by P and 2) that the local vapor pressure, P , can be represented 
n v, x 

by P . Subsequently, we shall see that P (^P ) cannot differ 

n, T n, T v, x 

appreciably from P at (or near) the leading edge of the annular cavity. 

Then the pressure difference, P^ ^ - P , is the driving potential for 

vaporization, and P - P causes condensation to occur. 

n n, T 

To graphically illustrate the foregoing discussion, we compare 

pertinent pressure profiles on figure 5. 8. This plot shows the venturi 

pressure distributions for non-cavitating flow (P ) and cavitating flow 

(P and P ) ; the venturi inlet liquid bulkstream vapor pressure, P , 
n n, T r v 

is also shown. All of these data are considered typical and were ob- 
tained from Run 123B, as tabulated in appendix A. The data plotted on 
figure 5. 1 can be obtained directly from the information shown on 
figure 5. 8. P , on figure 5. 8, was computed by combining the test 
data for Run 123B with the pressure coefficient data given on figure 3. 5. 
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lgure 5. 8 Typical pressure distribution through test section 
for cavitating and non-cavitating flow. 



P and P , on figure 5. 8, were plotted directly from test data for 
n n, T 

stations 2 through 9; upstream of station 2, P and P were estimated 

n n, T 

The shape of the P curve upstream of station 2 is in agreement with the 

generalized results of Rouse and McNown [25]; i. e. , for various cavi- 

tating bodies [25], as cavitation progresses, the minimum pressure 

increases while shifting slightly downstream and forming a shallower, 

but wider pressure 'trough'. The P and P curves are coincident just 

n x 

upstream of the minimum pressure point. The shape of the P^ ^ curve 

upstream of station 2 is compatible with visual and photographic observa 

tions from this study and previous work [15], These visual data indicate 

that the leading edge of the cavity always originates near, or slightly 

downstream of the non- cavitating minimum pressure point. If we select 

the cavitating minimum pressure point as the leading edge of the cavity, 

we may construct the P m (^P ) curve upstream of station 2 as 

follows: 1) the local vapor pressure, P , is coincident with P in 

n, T v 

the venturi inlet, 2) P increases slightly at the stagnation point on 

n, T 

the quarter-round contour due to the increase in P (or P ), 3) P m 

n x n, T 

decreases as P decreases to a minimal value and 4) P m decreases 
n n, T 

rapidly, as a result of vaporization, to the measured value of P . 

2, T 

The shaded regions between the P and P curves on figure 5. 8 indi- 

n n, T 

cate where apparent vaporization and condensation occur. The shaded 
region upstream of station 2 indicates that P - P should provide ade- 
quate estimates for cavity growth in this area. From figure 5. 7, we 
recall that apparent thermodynamic equilibrium generally exists at 
stations 2 and 4. Then little or no growth would be expected in this 
vicinity, see figure 5. 8; however, from physical considerations, the 
cavity thickness must increase between stations 2 and 4 in the venturi 
throat. All cavitating bodies [25] show cavity thickness increasing 
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with increasing cavity length in the front portion of the cavity. Although 
the instrumentation is not sufficiently precise to indicate the exact mag- 
nitude of P - P , it is apparent that some additional vaporization 
n, T n 

must occur in this region to support the increased cavity thickness. 

The measured values of P and P as shown on figures 5. 1 to 5. 5, 

n n, T 

5. 7, and 5. 8, imply that only slight vaporization can occur between 
stations 2 and 5. The above mentioned figures also indicate that transi- 
tion from vaporization to condensation generally occurs between stations 
4 and 6 and commonly occurs between stations 4 and 5. Then, from this 
discussion we may conclude that vaporization normally prevails upstream 
of station 5 and condensation occurs downstream of this location. 

It is interesting to note that station 5 coincides closely with the 
venturi throat-to-diffuser transition, see figure 5. 8. Since the diffuser 
was designed to provide efficient liquid pressure recovery, a condition 
favoring condensation, it is reasonable that condensation should com- 
mence near station 5. It appears that our data, e. g. , see figure 5. 8, 
are in good agreement with physical considerations. 

To avoid confusion, we emphasize that: 1) P m P ) differs 

n, T v, x 

appreciably from P^ at all cavity locations downstream of station 2, 

see figure 5. 8, 2) use of P - P in eq (5. 2-1) is valid only at locations 

v n 

upstream of station 2 and therefore 3) positive values of P - P down- 

v n 

stream of station 2 do not predict that vaporization should extend over 

the entire length of the cavity. Likewise, it should not be erroneously 

concluded from inspection of eq (5. 2-2) that the cavity will not collapse 

unless P exceeds P _ . Measured values of P and P m were chosen 
n n, T n n, T 

to approximate P^ ^ and P^ respectively. The actual position of the 

P p and P profiles on figure 5. 8 are unknown; therefore the unknown 
%x v,x 

quantity P. - P can dictate cavity collapse, even though the cavity 
^ , X V , x 

vapor is in apparent thermodynamic equilibrium, i. e. , P = P . 

n n, T 
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In summary, formulae applicable to spherical bubbles are being 
used to estimate vaporization and condensation rates at specific loca- 
tions within an annular cavity; the bubble formulae apply to a uniform 
pressure and temperature field (liquid), but the annular cavity is im- 
mersed in a nonuniform pressure and temperature field. Pressure and 
temperature data from this study are used in conjunction with these 
formulae to perform computations at fixed points within the annular 
cavity. In this way, the simplified formulae are adapted to our com- 
plex problem (where heat transfer and pressure variations prevail). 

The cavity growth rate is bounded by eq (5. 2-1), but the collapse 
rate increases without bound as the cavity collapse progresses. Care- 
ful study of eq (5. 2-2) indicates that R- increases almost linearly with 

VJ 

decreasing R until R ^R./2; thereafter, R^ increases according to 

-3/2 1 ^ 

the R relationship. Therefore, it can be shown that the additional 

time required for completion of cavity collapse after R = R. /4 is neg- 
ligible. More than 97 percent of the total collapse time is incurred in 
reducing the bubble radius to R. /4. We may then obtain a character- 
istic collapse rate by integrating eq (5. 2-2) over the radius interval, 
R./4 to R, , and calculating the mean collapse rate, R : 

li ^ 



(5. 2-3) 


The ratio of apparent collapse/growth rates may then be written as 


0. 5 


V r g 


2. 5 


n n, T 

P rr, - P 
n, T n 


(5. 2-4) 
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So, for equal vapor and liquid metastabilities, i. e. , ^ 

= P - P , a bubble will collapse 2. 5 times as fast as it will grow, 
n, T n 

In our experiments, P^ ^ - P^ at the leading edge of the cavity far 

exceeds P - P in the middle and aft regions of the cavity. Simulta- 

n n, T 

• • 

neous inspection of eq (5. 2-4) and figure 5. 8 indicates that R /R can 

C G 

vary widely with the values selected for P - P and P - P . 

n n, T n, T n 

For the purpose of this discussion, we simply need the average values 

of these pressure differences on either side of station 5. We can then 

compare, on the average, the growth rate upstream of station 5 to the 

collapse rate downstream of station 5. Referring to figure 5. 8, and 

using the data tabulated in appendix A, we obtain for these average 

2 

values, P _ - P 3. 8 p si (2. 6 N/cm ) and P - P m « 1. 5 psi 
n, T n n n, T 

(1. 04 N/cm?). Data uncertainty and data near the irregular, trailing 
edges of the cavity (station 9) are neglected in this estimate; actually, 
station 9 data may be included with little effect on the computations 

performed herein. Applying these mean values of pressure-difference 

• • 

to eq (5. 2-4), we obtain R /R 1. 6, i. e. , the growth time exceeds 

C G 

the collapse time by about 60 percent. Applying this same approach 

• « 

to figures 5. 2 through 5. 5, we find that R /R varies from about 

C G 

1. 5 to 2. 5. 


Within the bounds of this rough computation, we may then con- 
clude that the growth and collapse rates are not much different. For 
constant liquid velocity in the cavitated region, this implies that vapori- 
zation and condensation should occupy about the same fraction of cavity 
length. Referring to figure 5. 8, we see that the axial distance from 
the minimum pressure point to station 5 is about half the distance 
between stations 6 and 8 and about one-third the distance between 
stations 6 and 9, For shorter cavities, this distance ratio will approach 


unity and our computations will agree more closely with the test data. 
Because of the assumptions and approximations involved, the analysis 
presented here is highly simplified, but it certainly does not detract 
from the results shown on figure 5. 7. It seems totally plausible that 
metastable vapor should exist in the downstream regions of the cavity. 

To conclude this conjectural discussion, we must estimate the 
liquid velocity at the liquid- vapor interface. We will assume that the 
area available for liquid flow in the cavitated region is constant and is 
equivalent to the venturi throat area. For the experiments reported 
here, the liquid velocities in the throat varied from about 150 to 300 ft/s 
(45. 7 to 91. 5 m/s)--these velocities were undoubtedly higher due to 
the presence of the annular cavity of unknown thickness. The time 
required for a particle of liquid to traverse the liquid-vapor interface 
of the longest cavity is then computed by dividing the cavity length by 
the velocities given above; 0. 8 to 1. 6 milliseconds transit time is 
required. This transit time is shorter for higher liquid velocities and 
shorter cavities. Obviously, little time is available for vaporization 
and condensation to take place, and the metastable behavior depicted 
on figure 5. 7 is not too surprising. 

Actually, the thermal- response of the liquid, to rapidly varying 
pressure, as a particle of liquid traverses the venturi contour will 
determine where vaporization effectively stops and condensation starts. 
The transition from vaporization to condensation will occur at slightly 
different locations within the cavity as the fluid temperature, velocity, 
cavity length, etc. , are varied. An analytical model, with appropriate 
provisions for heat transfer and time -dependent pressure, is required 
to shed more light on this topic; at this point, we conclude our 
conjecture. 
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With respect to the foregoing analysis, we note that the liquid 
velocities in the vicinity of the cavities exceed the two-phase sonic 
velocity for hydrogen. We should also note that on occasion we have 
observed, on film, apparent n shock waves" traveling back and forth 
within the cavitated regions. These "shocks" manifest themselves as 
sharply- contrasted density discontinuities, made visible by the accom- 
panying variation in refractive index. 

5. 3 Data Correlation 

The correlative technique developed by Gelder, et al. [15], was 

used to correlate the data from this experiment. This technique has 

Z 

proven highly successful and is an extension of the B-factor concept. 

A brief history on the development of the B-factor concept and a simpli- 
fied and improved method for computing B-factor were recently pub- 
lished [20], Although the current B-factor approach is not entirely 
compatible with the physical processes of cavitation, it is well estab- 
lished, adequately documented, and provides good results [4, 5]. 

Recent efforts [23] to obtain a more compatible theory have not yet 
resulted in a predictive technique that is less complex or less dependent 
upon experimental data. The basic correlation expression is based 
upon dynamically and geometrically similar cavities. Thus, for con- 
venience, we will frequently refer to this correlative expression as 
the ’similarity 1 equation. The similarity equation is used to correlate 
developed-cavitation data in similar test items and to predict the cavi- 
tating performance of a test item from fluid-to-fluid and from one 
temperature to another, when limited data from a single fluid are 


B is defined as the ratio of the volume of vapor to the volume of liquid 
involved in sustaining a developed vaporous cavity. 
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available. The similarity equation in its final form is given [3] as 


B = (B) 


ref 


X V o,ref' V ref ' X o,ref' 


(5. 3-1) 


the symbols are identified in the nomenclature of this paper. To 
account for differences in theory and experiment, the exponents on the 
various terms in eq (5. 3 -1) are individually evaluated using the experi- 
mental data and theoretical data [20] for B as follows: 

(1) A theoretical value of B is obtained for each experimental 

data point using the measured cavity pressure-depression 

(P - P_), T , and the calculation method outlined in reference 
v 2 o 

[20]. This value of B is derived under the assumption that the 
vaporous cavity is formed by the isentropic vaporization of liquid 
and is referred to herein as B^ . 

(2) One experimental data point is arbitrarily chosen as a 

"reference"; the V D , and B from the chosen test are 

o o t 

then inserted into eq (5. 3-1) as constants where the subscript 
M ref M occurs. As explained in appendix C, a reference data 
point is ultimately chosen as that run which provides the best 
solution to eq (5. 3-1) and is therefore most representative of 
all data being correlated. 

(3) Values of ot 9 V 9 and D from each data point are then 

o o 

inserted into eq (5. 3-1) as the non-subscripted terms. This 
produces an equation for every data point except the one chosen 
as a reference. Note that the unknowns in the set of simultaneous 
equations are B and the exponents El, E2, E3, and E4. 

(4) The digital computer is then programmed- - see appendix C-- 
to determine the values of El, E2, E3, and E4 that provide the 
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best agreement between the computed B (left side of eq (5. 3-1)) 
and the theoretical obtained in step (1) for all data points con- 
sidered. The least squares technique used to evaluate these 
exponents is by no means trivial and is treated with appropriate 
detail in appendix C. This method ensures that the calculated 
B values and the B^ values for each data point are as nearly iden- 
tical as possible; the exponents computed in this manner represent 
the best possible agreement between experiment, eq (5. 3-1), and 
the isentropic flashing theory. Recall that the theoretical B^ 

(step 1) and the calculated B (eq (5. 3-1)) both rely on experi- 
mental data at each data point. 

The approach outlined above is applicable to any number of terms 
in eq (5. 3-1). As shown in the computer programs listed in appendix C, 
two additional terms were provided. Kinematic viscosity was included 
as a correlating parameter, because it has shown promise in some 
recent studies [23]; also, modification of the Gelder, et al. [15], theory 
to account for convection introduces the viscosity term. A surface ten- 
sion term was also tried, because this fluid property has long been rec- 
ognized [24] as a candidate correlating parameter and because the Weber 
number is vital to dynamic similitude studies involving the formation of 
bubbles, break-up of liquid jets, etc. Both of these fluid properties 
could justifiably be introduced through a dimensionless analysis approach; 
viscosity entering through the Prandtl or Reynolds numbers and surface 
tension by the Weber number. 

In eq (5. 3-1), the fluid physical properties are evaluated at P 
and T^ and the standard deviation in B is computed for each set of 
exponents; the individual exponents may be held constant or chosen by 
the computer. The standard deviation in B factor is minimized in the 
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computer programs when one or more of the exponents is selected by 
the computer; the absolute minimum standard deviation is obtained when 
all of the exponents are selected by the computer. The standard devia- 
tion is simply computed in those cases where the exponents are held 
constant. The set of exponents that produces minimum standard devia- 
tion in B is selected as the best correlative solution for any particular 
batch of data; i. e. , the standard deviation is a measure of the validity 
of the similarity and isentropic-flashing theories, as both are evaluated 
from experimental data. Data from the initial study [16] and the new 
data were separately and collectively correlated using the approach 
described herein. The results of this effort are given in table 5. 1, 
along with the results of others [3], and are discussed in the following 
section of this paper. 

5. 4 Discussion of Correlative Results and Data 

The similarity equation, used to correlate cavitation performance 
of a particular flow device from fluid -to -fluid, was fitted with numeri- 
cal exponents derived from the experimental data of this study. The 
similarity equation and exponent data for hydrogen- refrigerant 114 were 
obtained from the literature; the numerical exponents for hydrogen- 
refrigerant 114 are compared in table 5. 1 with those deduced from this 
experiment and the initial study [16]. The exponents given in table 5. 1 
were obtained with a least- squares fitting technique and a digital com- 
puter; the suitability of the various exponents to the experimental data 
is indicated by the standard deviation in B-factor as explained previously. 
In this study, the value of B ranges from two to five. 

As indicated in table 5. 1, there is little difference in the exponents 
obtained from the old data (line 1), new data (line 2), and combined old 
and new data (line 3); also, the data of Moore and Ruggeri [3] for 
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Table 5. 1: Correlation of liquid hydrogen data using the 'similarity' equation. 


I 
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points (including "ref" data point), B is computed from isentropic-flashing 


hydrogen- refrigerant 114 is in reasonable agreement, excluding the 
diffusivity term. The lack of variation in a ( < 10 percent) explains why 
El tends to a negative number when correlating with liquid hydrogen 
alone. Gelder, et al. [15], obtained a value of 0. 5 for El when cor- 
relating refrigerant 114 data with 25 percent variation in There 
was over 400 percent change in or i n the hydrogen-refrigerant 114 data 
correlated by Moore and Ruggeri [3], and thus the value for El reported 
in line 4 of table 5. 1 is to be preferred. The mathematical technique 
used to derive the exponents can easily pick an extraneous value for any 
of the exponents if there does not exist significant variation in the cor- 
responding physical parameter. 

Introduction of the viscosity and surface tension terms into 
eq (5. 3-1) did not significantly improve the hydrogen data correlation, 
and consequently exponents for these terms are not included in table 5. 1; 
however, these terms may be important correlating parameters when 
attempting to correlate data from fluid -to -fluid. It is anticipated that 
both of these parameters may improve data correlation where sufficient 
variation in the fluid properties occurs. In the hydrogen data of this 
study, the viscosity varied less than 20 percent and the surface tension 
varied less than 40 percent. Correlation of the hydrogen-refrigerant 
114 data would most likely be improved by using one or both of these 
terms. 

The diameter term in eq (5. 3-1) was not included in table 5. 1, 
because our tests were conducted with only one venturi size. Moore 
and Ruggeri [3] obtained an exponent value of -0. 1 for the diameter 
term, based on tests using refrigerant 114 in two different venturi 
sizes. Those tests were performed with a venturi identical to the one 
used in this study and with a larger (1. 414:1) geometrically similar 
venturi. 
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The arithmetic mean value of the developed cavitation parameter, 

K . , does not vary appreciably for all data considered in table 5. 1; 
c , min _ 

this is an important result, and constant K . , eq (5. 3-1), and the 

c, mm 

isentropic-flashing theory are used to predict [4, 5] the cavitating per- 
formance of equipment. It is anticipated that nature will triumph and 

K . for other cavitating models will not remain constant for all 
c, mm ~ __ 

fluids, cavity lengths, velocities, temperatures, etc. K . , of 

c , min 

course, varies widely with model or equipment geometry, as does the 
pressure coefficient, C^. Also, the conventional cavitation parameter 
for developed cavitation, K , varies with flow conditions for any particu- 
lar geometry, e. g. , see table A- la and Rouse and McNown [25]. 

The experimentally-determined exponents in eq (5. 3-l)--as listed 
in table 5. l--are based upon pressures measured within the cavity, 
near the leading edge where equilibrium prevailed in all fluids tested 
to date. It is clear that the experimental conditions under which the 
exponents were determined are compatible with the isentropic-flashing 
theory and local equilibrium throughout the cavitated region is not 
required. Also, it appears plausible that exponents derived in this 
manner could easily mask the effects of slight metastabilities (or per- 
haps unstable equilibria) within vaporous cavities. 

4 

The data given on figures 5. 1 to 5. 5 indicate that some of the 
cavities were shorter than their visual (as observed on film) length. 


K . for every hydrogen data point was within 7 percent of K 

c, mm c, mm 

(excluding Run No. 117). 

4 

The pressure-depression should be near zero at the trailing edge of 
the cavity and the actual length of the cavity may be estimated by extrapo- 
lating the data to zero pressure-depression. 
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The actual length of the cavity and the observed length differ because of 
irregular trailing edges of the cavity and the difficulty in judging the 
visual length. Both actual and visual cavity lengths were used to corre- 
late the data, and they produced essentially the same results. 

The pressure depression in the cavitated region is determined by 
subtracting the measured cavity pressure, in one case, and the satura- 
tion pressure associated with the measured cavity temperature, in the 
other case, from the vapor pressure of the liquid entering the test 
section. 

In the data reported here, the minimum measured cavity pressure 

was less than bulkstream vapor pressure by as much as 15. 13 psi 
2 

(10. 44 N/cm ); these pressure-depressions are obtained by subtracting 

from P in the tabulated data of appendix A. 

2 v 

Cavity pressure-depression increases with increasing cavity 

length, fluid temperature, and velocity for these tests. These trends 

are graphically demonstrated in figures 5. 1 to 5. 5 and in reference [16]. 

The data in appendix A also readily disclose these trends if we permit 

only one of the three parameters (£, V , and T ) to vary at a time and 

o o 

note the value of P - P^ for two different values of the parameter being 

varied. As an example, T^ and are relatively constant in Runs 125A 

and 125B, but the cavity length differs by a factor of three; P - P^ 

for the longer cavity (Run 125B) is about four times the P - P^ for 

v 2 

the shorter cavity (Run 125A). Similar comparisons may be made for 
variations in and T , or portions of the data may be plotted collec- 
tively on single graphs, etc. , to demonstrate these trends. Typical 
plots (14 graphs) are given in reference [16]. 

Because of its popularity in the pumping machinery field, pressure- 
head has been included in the data tabulated in appendix A; however, 
provision of these data does not mean that the authors approve of the 
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use of pressure -head terms. On the contrary, this practice is to 
be discouraged [20] in computations related to the cavitation process. 
Mathematical conversion of pressure to pressure-head merely requires 
evaluation of the liquid density at the point of measurement; however, 
selection of the appropriate liquid density can be a bit perplexing. Fig- 
ures 5, 1 to 5. 5 indicate that the measured pressures and temperatures 
at the same axial location, within the cavities, are generally not in 
stable thermodynamic equilibrium. Also, due to the thermal expansivity 
of liquid hydrogen, the bulkstream temperature changes appreciably as 
the liquid flows through the venturi. The following methods were used 

to calculate pressure-head from the cavity measurements: (1) Head (h ) 

n 

was calculated from measured cavity pressure by using the saturation 

density at the measured pressure. (Z) Head (h ) was calculated from 

n, T 

measured cavity temperature by using the saturation density at the 
measured temperature. Both values of head are given in the tabulated 
data in appendix A. 

The cavitation parameter for fully developed cavitation, K , was 
calculated and tabulated for each run, see appendix A. 
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6. CONCLUDING REMARKS 

Pressure and temperature profiles were measured within fully- 
developed, vaporous hydrogen cavities; the results are given in appen- 
dix A and on figure 5. 7. In general, the measured pressure and tem- 
perature depressions were not in stable thermodynamic equilibrium; 
the pressure-depressions obtained from the cavity temperature meas- 
urements are usually greater than those derived from the measured 
cavity pressures, see typical plots in figures 5. 1 - 5. 5. Specifically, 
the experiments indicate that the cavity vapor is normally in thermo- 
dynamic equilibrium near the leading edge of the cavity, while con- 
siderable thermodynamic metastability occurs in the central and trail- 
ing regions of the cavity. This behavior is attributed to lag in the 
thermal-response of the liquid, to rapidly varying pressure, as a par- 
ticle of liquid traverses the test section contour. Also, it is shown by 
simplified analysis that inertial effects alone could at least partially 
account for this behavior. This study confirms the metastability phe- 
nomena reported earlier [16]. 

The experimental data from this study and previous work [16], 
were used to fit a 'similarity 1 equation with numerical exponents, see 
table 5. 1; favorable comparisons are drawn between the old data [16], 
new data, and the work of others [3]. This 'similarity' equation, 
coupled with the isentropic flashing theory [20], has been proven useful 
in predicting [4, 5] the cavitating performance of liquid pumps from 
fluid-to-fluid. Introduction of viscosity and surface tension parameters 
into the 'similarity' equation did not improve the correlation for the 
hydrogen data; however, it is anticipated that both of these parameters 
may improve fluid-to-fluid correlations. An attempt should be made 
to update the refrigerant 114 -hydrogen correlation [3] using these two 
additional parameters. 
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Recent advances in low temperature thermocouple thermometry- 
permit more precise, local temperature measurements within vapor- 
filled cavities; consequently, experimental data are sufficiently reliable 
to evaluate future modifications in cavitation modeling theory. 

7. ^NOMENCLATURE 

B = ratio of vapor to liquid volume associated with the 

sustenance of a fixed cavity in a liquid 

BCONV = symbol used in computer programs for B (left side 

of eq (5. 3- 1 ) ) 


B FLASH 


B 

t 

C 

P 


c 

p 

D 

o 


symbol used in computer programs for B^_ 

B derived from isentropic flashing theory (Ref. [20]) 

2 

pressure coefficient [ = (h - h Y/(V /2g )1 

x o^o c 

v 2 

minimum pressure coefficient [ =(h - h Q )/(V 
test section inlet diameter 

conversion factor in Newton’s law of motion (gravita- 
tional acceleration) 


h 

n 


(n = 2, 4, 5, 6, 7, 8, or 9): head corresponding to cavity 
pressure, measured at a particular instrument port 
in wall of plastic venturi 
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h 

O 


h 

v 


h 

X 


K 

c, mm 
K 

c, mxn 


K 


v 


l 

P 




X 


(n = 2, 4, 5, 6, 7, 8, or 9): head corresponding to the 
saturation pressure at the cavity temperature, meas- 
ured at a particular instrument port in wall of plastic 
venturi 

test section inlet head corresponding to absolute inlet 
pressure 

head corresponding to saturation or vapor pressure at 
the test section inlet temperature 

head corresponding to absolute pressure, measured 
at wall of plastic venturi at distance x, downstream 
of the minimum pressure point--for non- cavitating 
flow 

head corresponding to the minimum absolute pressure 

on quarter-round contour of plastic venturi, computed 

v 

from expression for C 

P 

developed cavitation parameter, based on minimum 

2 

measured cavity pressure f = (P - P )/(P V /2g )] 

o Zoo c 

arithmetic mean value of K . for a set of data 

c, mm 

points 

developed cavitation parameter [= (h^ - h )/(V /2g )] 

visual (observed on film) length of developed cavities 

local bulkstream liquid pressure at distance x along 
the cavity interface -- variable for the annular cavity, 
but taken as constant in the development of eq ? s 
(5. 2-1) and (5. 2-2) 
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(n = 2, 4, 5, 6, 7, 8, or 9): absolute cavity pressure, 

measured at a particular station or instrument port 
in wall of plastic venturi- - also used to approximate 
in section 5.2 

ts, x 

(n = 2, 4, 5, 6, 7, 8, or 9): saturation pressure corre- 

sponding to the measured cavity temperature at a 
particular station or instrument port in wall of plastic 
venturi-- also used to approximate in section 5. 2 

test section absolute inlet pressure 

saturation or vapor pressure at test section inlet 
temperature 

local vapor pressure of liquid at distance x along the 
cavity interface- -variable for the annular cavity, but 
taken as constant in the development of eq's (5. 2-1) 
and (5. 2-2) 

absolute pressure, measured at wall of plastic venturi 
at distance x, downstream of the minimum pressure 
point- -for non-cavitating flow 

bubble radius 

differentiation of R with respect to time 

Reynolds number, based on test section inlet diameter 

’initial 1 bubble radius at the onset of collapse 

(n = 2, 4, 5, 6, 7, 8, or 9): measured cavity temperature 
at a particular station or instrument port in wall of 
plastic venturi 



T 

o 


V 

o 

X 


Greek 

Oi 

P 

o 

Subscripts 

C 

G 

ref 


bulkstream temperature in degrees Rank in e (Kelvin), 
of liquid entering the test section 

velocity of test liquid at inlet to test section 

distance measured from minimum pressure point on 
quarter-round contour along axis of plastic venturi 


thermal diffusivity of liquid at inlet to test section 
absolute viscosity of liquid at inlet to test section 
density of liquid at inlet to test section 


denotes state of bubble collapse 
denotes state of bubble growth 

reference run (data point), or test conditions, to 
which a computation is being referenced when 
attempting to correlate cavitation performance 
via the similarity equation (5. 3-1) 


Superscripts 


El 

E2 

E3 

E4 


exponent on thermal diffusivity ratio in eq (5. 3-1) 

exponent on test section inlet velocity ratio in 
equation (5. 3-1) 

exponent on cavity length ratio in eq (5. 3-1) 


EN 
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exponent on test section inlet diameter ratio in 
eq (5. 3-1) 

any one of the exponents used in eq (5. 3-1) 
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Table A-la. Experimental developed-cavitation data in venturi using liquid hydrogen (English Units) 


Appendix A: Experimental developed-cavitation data in venturi using liquid hydrogen. 
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Table A- la. (Continued) 
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Table A- la. (Continued) 
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Table A- lb. Experimental developed - cavitation data in venturi using liquid hydrogen (SI Units) 
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Table A- lb. (Continued) 
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Table A- lb. (Continued) 


I 
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Appendix B: Thermocouple fabrication procedure 

I. Support tubes 

A. Use 1 /16-inch dia. x 0. 012-inch wall stainless steel tubing. 

B. Cut 1. 25 -inch lengths of this tubing with jewelers file. 

C. Hold a tubing length in pin vise and deburr tubing with 

1 /16-inch diameter drill held in another pin vise, 

D. Soak these support tubes in trichloroethylene to clean. 

II. Thermocouple junctions 

A. Skin thermocouple wire insulation with #320 abrasive cloth. 

5 

B. Tin the Chromel wire using stainless steel soldering flux 
and solid core (50/50) solder. Wash off flux. 

C. Twist Chromel-Au (0. 07 at. percent Fe) wires together to 
form junction and solder using resin-core radio solder. 

Use voltage-controlled soldering iron and do not get Au 
wire too hot--use just enough heat to melt solder. 

D. Clean soldered junctions with a solvent, and then coat all 

5 

uninsulated surfaces with G. E. Varnish. 

III. Support tube and thermocouple junction assembly 

A. Wipe wires with solvent (avoid varnish). 

B. Thread junctions and suitable length of extension wire into 
1 /16-inch diameter support tubes. 

C. Attach a transparent vacuum hose to the junction (sensor) 
end of the assembly. 

5 

D. Mix epoxy: Emerson and Cuming, Inc. , Stycast 2850 GT 
with room temperature cure activator. 


Identification of a manufacturer and a manufacturer’s product has 
been necessary to make the results of this work meaningful and in no 
way implies a recommendation or endorsement by the National Bureau 
of Standards, or by the National Aeronautics and Space Administration. 
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E. Apply epoxy to the non- junction end of the support tube; 
then use the vacuum-assist and manipulate the thermo- 
couple wires in the tube until epoxy appears at the junction 
end (inside transparent vacuum hose) of the support tube* 

F. Remove vacuum hose and excess epoxy from both ends of 
the support tube. 

G. Make final placement of junctions and allow adequate cure 
time. Junctions should project slightly from the support 
tubes and should be visibly centered in the support tube; 

i. e. , the thermocouple wire should not be in contact with 
the support tube on either end. 

IV. Installation 

A. The 1/16 -inch diameter thermocouple probes are then 
inserted into the plastic venturi through compression 
fittings as indicated in figure 3. 6. 

The reference thermocouples (figure 3. 7) were fabricated in the 
same fashion, except that all of the 1/16-inch diameter support tubes 
were first soldered, with parallel axes, inside of a single 1 /4-inch 
diameter support tube. This precaution assures that all of the ref- 
erence thermocouple junctions will be at the same temperature. 



Appendix C: Computer programs for correlation of 

developed cavity data 

This appendix outlines the details of the two digital computer 
programs used to determine the exponents on eq (5. 3-1). The pro- 
grams will be described, in turn, and will be referred to as (1) the 
variation-of-parameter solution and (2) the Bjorck solution. The Bjorck 
solution is a sophisticated n least-squares M method, and the variation- 
of-parameter solution is a basic least- squares approach. The latter 
program uses the method of "steepest descent" to locate a minima in 
the selected ordinate--in this case, standard deviation in B-factor. 

The Bjorck program always finds the absolute minimum ordinate, while 
the variation-of-parameter program may find a local minima; however, 
the method of steepest descent tends to locate the absolute minimum. 

The Bjorck program consumes little machine -time and simultaneously 
determines, without bound, all of the exponents in eq (5. 3-1). The 
parametric program is much slower and incrementally varies one 
exponent at a time, conveniently providing realistic bounds on the ex- 
ponents. Providing bounds for the exponents is a necessary feature 
when correlating data with only a slight variation in one or more of the 
physical parameters, e. g. , <*, V^, or & m 

To illustrate this point, consider a group of data where a varies 
by only 5 percent and the data scatter is about 5 percent--the Bjorck 
program may pick El = 1. 0 or -1. 0, but by selecting the incremental 
range and initial value of El, we can force the variation-of-parameter 
program to pick a value for El between any limits that we choose, e. g. , 
between 0.2 and -0.2. Because the variation-of-parameter solution 
may not find the absolute minimum ordinate, the Bjorck program pro- 
vides a ready check on the parametric solution. 
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For data with appreciable variation of all physical parameters, 
the Bjorck program is far superior, because it is much faster and more 
reliable; therefore, this program is ideal for most fluid -to -fluid cor- 
relations. The variation-of-parameter solution is usually best for 
single fluid correlations where the physical parameter data may not 
vary much. 

To feel confident that the best solution has been obtained, one 
must possess a basic understanding of the problem and exercise con- 
siderable intuition, i. e. , curve-fitting is still an artj While the digital 
computer is a powerful tool that can rapidly solve complex and cumber- 
some problems, it is restrained by the mathematical techniques at our 
disposal. The method of 'least- square s' is no exception, and one should 
be aware of certain pitfalls ; some of the hazards of attempting to blindly 
fit data with this technique will be revealed herein. There are two 
practical reasons why the least squares method [26] is used extensively 
for data-fitting problems: (1) the condition of minimizing the sum of 

the squares of the deviations is a convenient and efficient method of 
handling large quantities of data with computer routines, and (2) sta- 
tistical theory claims that when certain necessary least squares con- 
ditions are met, the least squares solution is the 'best' obtainable. 
Variation-of-parameter solution : This program uses the least squares 

technique and the method of steepest descent to locate a minima. To 
fully understand how the computer program is formulated, we rewrite 
eq (5. 3-1) in abbreviated form, 

B = B ref (°ref /a ) ( W ref ) ’ ' ' > (C - 1) 


61 


L 



Note that in the computer program, see page Cll, B = BCONV, 

B = B = BFLASH, for the selected reference run, ot = ALPHA 
ref t ref 

(REF), etc. In the discussion that follows, computer-program symbols 
will always appear in capital letters and the symbol, *, will denote 
multiplication. The method of least squares is applicable to non-linear 
problems if they can be converted into equations where the parameters 
enter linearly. Equation (C-l) can then be handled by this technique if 
it is linearized and if we can reasonably estimate an approximate solu- 
tion, El = E1,0, E2 = E2,0, etc. Wylie [26] shows that the deviation 
in B is given by a Taylor's series expansion about the approximate solu- 
tion point. 


6 B 


af 

3E1 


❖ 6 El + 
El, 0 


af 

dE2 


* 6 E2 + . 
E2, 0 


etc. , 


(C-2) 


where B = f (El , E2, E3, . . . , etc. ) and 6 E 1 to 6 EN is the exponent 
increment used in the computer program. Recall that B^ , a > V , 
etc. , are input data and the exponents, El to EN are to be determined. 

Performing the indicated differentiation on the first two terms 
of (C-l), we obtain 


6B 


El 


E2 


B 


ref 


+ 


(¥) (^) •■(¥) 

x o, ref * 


6E2 . 


This expression then reduces to 

5 B = BCONV | ln[ — ] *6El+£] 


■n *6El+£n(— ^ )* 

U ' V V o,ref/ 


6 E2 


, (C-3) 
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for the two terms considered. Now if we only permit one exponent to 
vary at a time, e. g. , El, we may write 


6 B 


BCONV 



6E1 . 


(C-4) 


Similar expressions may be written for each exponent (EN) and 
its accompanying log term. Using (C-4), we obtain an equation of con- 
dition for each data point, and every equation in this system is linear 
B f I 

in 6EN; the | q becomes the coefficients in this linear sys- 

tem of equations. The method of least squares is then applied directly 
to this system- -each equation of condition is multiplied by the coeffi- 
cient of the variable in that equation and summed [26] --to obtain 



Sf 

a en 



(C-5) 


This expression derives directly from application of the least squares 
method to eq (C-2), recalling that only one exponent is being varied. 

Upon completion of this computation, the estimated solution is corrected 
by an appropriate incremental value and the entire process cycled until 
the calculated and estimated exponent values converge. In this approach, 
we are minimizing on the expression in (C-5), a necessary condition for 
minimization of the variance in B. We may, in fact, combine (C-2) and 
(C-5) under the constraint of single exponent variation, to obtain 


E (6B)2 




(C-6) 


Equation (C-6) will be recognized as the familiar "propagation of 
variance" rule, and (C-6) is minimized in the conventional least squares 
solution; i. e. , variance in B is minimized. The term BCONV * £ n 
(Ratio) in (C-3) and (C-4) may then be viewed as an 'apparent weight 
factor', applicable to this particular problem where each exponent is 
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altered sequentially. For a numerical solution, we will increment one 
of the exponents at a time, until all of the exponents have been changed-- 
this process will be repeated several times to provide sufficient range 
for each exponent- -then this incrementing process is repeated in ever- 
decreasing incremental steps. Repeating this process a sufficient num- 
ber of times with appropriate sized increments results in the Tf best M 
solution possible. This program uses two criteria to evaluate the qual- 
ity of the data-fit, (1) MNBDIFPC = the mean percent difference in B 
and (2) STDDEV = the standard deviation in B. These two terms are 
defined as follows: 


MNBDIFPC = 


y 

j BDIF | * 100 ( 



d BFLASH ( 



(NPTS-1) 


STDDEV = 


Z(BDIF) 2 

(NPTS-1) 


1/2 


where 6 B = BDIF = BCONV-BFEASH and NPTS = number of data points. 
This variation-of-parameter program functions as follows: 

(1) All input data, including BFEASH, are read and stored in com- 
puter memory--statement 1. 

(2) The reference run is indicated, and initial estimates for each 
exponent are read and stored--whether the exponent is to be 
fixed or varied is determined-statement 3. 

(3) The physical parameter ratios are computed and stored-- 
statement 9. 

(4) These ratios are then exponentiated- -statement 23. 

(5) Each exponent is then sequentially incremented in ever-decreasing 
steps- -statement 11 + 1 . 

(6) STDDEV and MNBDIFPC are evaluated--statements 20 + 2 and 
27 + 2, respectively. 
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(7) Pertinent dimensionless numbers may be computed upon demand- - 
statement 30. 

(8) The average value of K .is computed- -statement 35 + 1. 

c , mm 

Pertinent information is printed out as it is generated during the 
program operation. An option is included that permits calculation of 
the arithmetic mean n reference data" --statement 51. This calculated 
''reference run" always provides good results and thus permits a com- 
parative evaluation with actual runs when they are used as the reference 
dataj 

Bjorck solution; This program, though sophisticated, is programmed 
to provide a conventional least-squares solution to data-fitting problems. 
Because this program simultaneously determines the exponents in (C-l), 
we must rewrite (C-l) in the conventional 'linear 1 form; taking the loga- 
rithm of both sides of (C-l) yields 


a 


f(B) = In 


= El*tn 


ref , 


ref 
O' 


V 


+ E2*tn 


V 


|+ . • • etc. 


o, ref 


Recall that the ’propagation of variance 1 rule requires that 

2 


£< 6£ > 2 =2:(!i * 6b ) 


(C-7) 


(C-8) 


Now that our original problem is written in ’linear 1 form (with 
logarithmic coefficients), minimization of (C-8) results in the ’best 1 
solution of (C-7) in its logarithmic or nonlinear sense. This solution 


requires that 


E( 


SB * 6 B ] be minimized, and thus 
is minimal. The latter term, though different, contains the essential 


ingredients of MNBDIFPC. To obtain the 'best 1 solution in the linear 
sense, i. e. , for minimum ^(6B) , we must multiply the left side of 
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leaving J^(6B) on the right side of (C-8). The 


(C-8) by 1 

2 * x 

(6B) is tlien minimized in the Bjorck program, STDDEV is minimal, 
and the result is the 'best' available in the linear sense of the problem. 

/ 3f \ 2 2 

Differentiating (C-7), we obtain l gg ] = 


1/B . This term is commonly 

called the weight function and ideally should be evaluated with BCONV; 

however, BCONV has not yet been computed, so we must use the isen- 

tropic flashing value, BFLASH, for the first computation. Thereafter, 

the weight factor could be improved by using BCONV; this nicety did not 

appear to be warranted and is not incorporated in the Bjorck program. 

Theoretically, the solution to our problem that produces the minimum 

2 

STDDEV is what we seek, and so the 'weight 1 function, (BFLASH) , 


should be used with all data. When this weight term is neglected, the 
Bjorck program will return results based on minimum ^ 
usually, these results do not differ much from those obtained for mini- 
mal STDDEV. The Bjorck program functions as follows: 

(1) All input data, including BFLASH, are read and stored in com- 
puter memory--statement 1. 

(2) The number of terms to be evaluated, the number of constraints 
on the solution, the weight function, and the reference run are 
all predetermined--statement 13. 

(3) The solution may be constrained so that any number of the 
exponents may be fixed- -statement 45. 

(4) The log ratios are computed and stored and the Bjorck sub- 
routine is called- -statement 66 for 'weight' =1.0 and statement 

2 

76 with weight = (BFLASH) 

(5) The exponents are returned from the Bjorck subroutine, and then 
BCONV, BDIF, STDDEV, MNBDIFPC, etc. , are calculated- 
statements 16 to Z0. 
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(6) Pertinent dimensionless numbers may be computed upon demand- - 

statement 32. 

(7) The average value of K .is computed- -statement 35 + 1. 

c , mm 

Pertinent information is printed out as it is generated, and again 
an option is included that permits calculation of "reference data"-- 
statement 81. These reference conditions are derived by calling upon 
the Bjorck subroutine to compute a representative value of BFLASH; 
this value of BFLASH is derived by first evaluating the basic B-factor 
equation to obtain values for El, E2, E3, etc. , for a set of data. The 
basic B-factor equation is simply (C-l) without the parameters sub- 
scripted with "ref. " Then the mean values of the physical parameters, 

01 , V , etc. , are computed. These mean values and the exponents, EN, 
o 

are then used in the basic B-factor equation to calculate a 'representative 1 
BFLASH- -statement 85 + 8. This 'representative' BFLASH and the mean 
physical parameters are then used as "reference data" for the Bjorck 
program. This calculation method yielded low STDDEY's for compari- 
son with those attainable using actual run data as the "reference data. " 
Hazards of using the Variation-of-Parameter program : Judicious selec- 
tion of the initial exponent values, though not vital, helps the computer 
program to converge on the best solution more rapidly. The initial 
exponent estimates should match the anticipated results as closely as 
possible. Care must be taken to choose appropriate incremental steps 
for the exponents; if the steps are too small, a local minima may be 
found, and if the steps are too large, the absolute minimum may be 
bypassed in favor of some other minima. The number of incremental 
cycles should be kept to a minimum to conserve computer time and still 
provide appropriate bounds on the exponents. Where the physical param- 
eter, e. g. ot 9 does not vary much, it is wise to check for exponent con- 
vergence; i. e. , be certain the exponent value is not continuously 
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'drifting 1 when the incrementing is concluded. Where the physical 
parameter does not vary much, the order in which the exponent is 
altered can drastically affect its final value, e. g. , if & varies less 
than ten percent and El is the first exponent altered, weird results 
may be returned. Causing El to be the second or third exponent to 
be altered will normally produce much more realistic values, i. e. , 
the exponent on a weakly varying parameter should not be one of the 
first to be altered. Another factor that can cause considerable con- 
fusion is the occurrence of the same physical parameter in more than 
one exponentiated term. As an example, it is desirable to evaluate & 

and kinematic viscosity, Mj/ P , in the same expression. Fluid density, 

o 

P , is inherent in both of these terms and can interact in such a way 
o 

that unpredictable results occur when both terms are used. Again, 
the order in which these terms are processed in the computer program 
may affect the outcome. Also, the initial selection of exponents can 
have some influence on this interaction as well as determination of the 
best solution. 

In this program and the Bjorck program, it is wise to have some 
means of hand- checking the results. A crude method of doing this for 
our problem consists of plotting graphs as follows: 

(1) Data points are selected so that only one physical parameter 
varies appreciably; 

(2) this parameter is then plotted as the abscissa, with B/B ^ 
plotted on the ordinate of log - log graph paper; 

(3) the best straight line is fitted to these plotted data; and then 

(4) the slope of this line is a rough estimate of the exponent value 
for that particular physical parameter, see eq (C-l). 
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In spite of all of the drawbacks of the variation-of-par ameter 
program, it is very useful and usually provides the lowest STDDEV 
and MNBDIFPC for a set of data points. Its major drawback is the 
computer time required, and its primary asset is the ability to pro- 
vide realistic bounds for the exponents. 

Hazards of using the Bjorck program; The major problem in using 
this program occurs where there is little variation in one or more of 
the physical parameters. In its present form, there are no provisions 
for setting realistic bounds on the exponents. The appropriate weight 
function must be used if minimization in the linear sense of the problem 
is desired. It is usually wise to run the program with and without the 
weight function, because minimization in the logarithmic or non-linear 
sense can provide superior results, i. e. , a lower value of MNBDIFPC 
and about the same STDDEV. This program requires little machine 
time and where the physical parameters vary significantly provides 
essentially the same results as the slower variation-of-par ameter 
program. It serves as a good check on the slower program, because 
within its mathematical limitations, the absolute minimum solution is 
guaranteed. Because of these niceties, this program is extremely 
useful in correlating data from fluid -to -fluid. 
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Variation-of-Parameter Computer Listing 


• JOB *2750453 »H0RD, 4 
• FTN»L,X,R 

PROGRAM BCONVECT 

DIMENSION DEL ( 3 ) .RUN (405) .DPM(405) ,DLIQ(405) »V0(405) .ALPHA (405) . 

1LENGTH (405) »DTS(405) .ROV I S ( 405 ) » ROS I G ( 405 ) » BFLASH ( 405 ) » NN < 405 ) . 

2EC5 > .RAT 1 0 ( 405 » 5 ) » RAT I OEXP ( 40 5 . 5 ) .BCONV(405 ) »BDIF<405) .HOLD (5 ) 

INTEGER HOLD, REF 

REAL KCMIN, LENGTH, MNBDIFPC 

DATA (DEL=0. 2,0. 02,0. 002) 

C THIS PROGRAM CORRELATES B-FLASH AND BCONVEC THEORIES — SIMILARITY 

C EQUATION — WITH EXPERIMENTAL DATA. 

C PROGRAM LOGIC9PROGRAM, DATA, BLANK CARD » EXPONENT-VALUE CARDS ( OPT I ONS 

C ) .BLANK CARD, MORE DATA ( E . G. D I FFERENT CAVITY LENGTHS AND PRESSURES- 

C — NN= 1 TO 6), BLANK CARD, MORE EXPONENT-VALUE CARDS, THEN EOF. 

C AFTER FIRST EOF ADD AN IBRANCH CARD THAT BRANCHES INTO DIMENSION- 

C LESS NO. CALCULATION, THEN ANOTHER EOF. THE REF RUN MUST BE SPECIFIED ON 

C THE EXPONENT CARD. IF THE 4TH OR 5TH EXPONENT 

C HAS INITIAL VALUE = 0 . 0 , THAT TERM WILL BE HELD ATI. 0, NOT USED, IN THE 

C SEEK-FIT OF THE OTHER EXPONENTS. THIS PROGRAM DOES FORCE-FIT INDIVIDUAL 

C TERMS .IF TERMS 1 AND 3 ARE TO BE HELD CONSTANT — TYPE 1 IN HOLD1 AND 

C 3 IN HOLD3, ETC. .LEAVE hold columns blank UNLESS A NON-ZERO exponent 

C IS TO BE HELD CONSTANT . HOLD TERMS APPEAR ON EXPONENT CARD. 

C DECK STRUCTURE . DATA, BLANK CARD , EXPONENT CARDS, BLANK CARD, DATA, 

C BLANK CARD, EXPONENT CARDS , EOF , I BRANCH CARD , EOF , EOF. 

C TO OBTAIN AN IBRANCH AFTER EACH BATCH OF DATA DATA, BLANK CARD, 

C EXPONENT CARDS, EOF, IBRANCH CARD, MORE DATA, ETC., 

41 1 = 1 

1 READ 1 0 1 , RUN ( I ) ,DPM< I ) , DL I Q ( I ) , VO ( I ) .ALPHA ( I ) , LENGTH ( I ) »DTS< I ) , 

1ROVISI I ) , ROS I G ( I ) , BFLASH ( I ) , NN ( I ) 

C A BLANK CARD IS REQUIRED TO END DATA. 

I F ( EOF ,60)40,42 

42 I F ( VO ( I ) » EQ. 0 .0 ) GO To 2 
1 = 1 + 1 
GO TO 1 

2 NPTS= I - 1 

THE FOLLOWING DO LOOP IS ELIMINATED WHEN ROSIG IS TO BE USED 

THIS DO LOOP IS NECESSARY TO EVALUATE SCALE EFFECTS WEBER NO. 

WILL BE MEANINGLESS WHEN THIS DO LOOP IS USED 
DO 44 1 = 1, NPT S 

ROS I G ( I ) = DTS ( I ) 

44 CONTINUE 

3 READ 102, ( E( I ) , 1 = 1,5 ) » (HOLDt I ) , 1 = 1 ,5 ) ,REF 
IF REF =1 THE FIRST CARD READ WILL BE USED AS REF RUN, IF REF=10 THE 10TH 
CARD WILL BE USED, IF REF=0 THE REF RUN WILL BE CALCULATED AND REF WILL 
BE SET EQUAL TO NPTS+1. 

I F ( EOF ♦ 60 >30,21 

21 I F ( E( 1 ) .EQ.0.0 ) GO TO 41 
IFfREF.EQ. 0.0)51, 24 
24 I F ( E ( 4 ) .EQ. 0.0)6, 5 

4 L I M=3 
GO TO 9 

5 I F ( E( 5 ) .EQ.0.0 ) GO TO 7 
GO TO 8 

C 6+1 STARTS CALC FOR FIRST 3 AND 5TH TERMS 
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6 I F t E C 5 ) .EQ.O.O) GO TO 4 
L I M = 5 

DO 22 1=1 »NPTS 

RAT 10 ( 1*1) = ALPHA t REF ) / ALPHA t I ) 

RAT 10 ( I *2 )=V0< I ) /VO ( REF) 

RAT 10 ( 1*3) =LENGTH ( I ) /LENGTH t REF ) 

RATI 0 ( I » 4 ) = 1 • 0 

RAT IOt 1*5) =ROS I G t I ) /ROSIG(REF) 

22 CONTINUE 
GO TO 23 

7 L IM=4 
GO TO 9 

8 LIM=5 
GO TO 9 

C STATEMENT 9 STARTS CALCULATION FOR FIRST 3»4, OR ALL 5 TERMS OF EQUATION. 

9 DO 10 1 = 1* NPTS 

RAT IOt I ,1)=ALPHA(REF)/ALPHA( I ) 

RAT IOt 1*2) =V0 ( I ) /VO (REF) 

RAT IOt I *3 ) =LENGTH ( I ) /LENGTH ( REF ) 

RAT IOt 1,4) =ROVI S ( I > /ROV I S ( REF ) 

RAT IOt 1 ,5 )=ROSIG( I ) /ROSIGt REF) 

10 CONTINUE 

23 DO 11 1 = 1, NPTS 

DO 11 J= 1 » L I M 

RAT IOEXP ( I , J ) =RAT I 0 ( I ,J)**E( J> 

11 CONTINUE 
M J=L I M 

DO 18 MM= 1 * 3 

DO 18 K= 1 > 20 

DO 18 J= 1 , M J 

I F ( HOLD ( J ) • EQ • J ) GO TO 18 

SUMBD I F=0 . 0 

DO 13 I = 1 * NPTS 

BCONVt I ) =BFLASH ( REF ) 

DO 12 J J = 1 » L IM 

BCONVt I ) =BCONV ( I ) *RAT IOEXP ( I »JJ ) 

12 CONTINUE 

BDIFt I ) =BCONV< I ) -BFLASH ( I ) 

SUMBD IF = SUMBDIF+BDI Ft I >*BDIF( I ) 

13 CONTINUE 

WRITE (61 >103) (E(I) *1=1*5) »SUMBDIF 
RES I DUE=0 .0 
DO 14 1=1 *NPTS 

EE=RAT IOt I*J) 

RESIDUE=RESIDUE+BDIF ( I ) *LOGF ( EE ) *BCONV ( I ) 

14 CONTINUE 
IF(RESIDUE) 15*18*16 

15 E(J)=E(J)+DEL(MM) 

GO TO 17 

16 E(J)=E( J)-DEL(MM) 

17 CONTINUE 

DO 18 1=1, NPTS 

RAT IOEXP ( I,J)=RATIO( I , J)**E< J) 

18 CONTINUE 


71 08 05 


71 



WR I TE ( 61 » 104 ) 

SUMBDIF=0.0 
SMBDIFPC=0.0 
DO 20 1 = 1, NPTS 

BCONV ( I )=BFLASH(REF) 

DO 19 JJ=1»LIM 

BCONV ( I >=BCONV( I ) *RATIOEXP( I »JJ) 

19 CONTINUE 

BDI F ( I ) =BCONV ( I ) — BFLASH { I) 

BDIFPC=(BDIF(I)/BFLASH( I) ) *100.0 
SMBDIFPC=SMBDI FPC+ABSF ( BDI FPC ) 

SUMBDIF=SUMBDIF+BDIF( I ) *BD I F ( I ) 

WRITE (61 >105 >RUN( I > , ALPHA < I ) * VO ( I ) , LENGTH ( I ) *ROVIS( I ) »ROSlG( I ) * 
1BC0NV ( I ) , BFLASH ( I ) »BDIF( I ) ,BD I FPC * SUMBD I F , DTS ( I ) *NN ( I ) 

C BDI FPC = PERCENT AGE BDIF AND SUMBD I F = SUMMAT I ON OF BDIF SQUARED 

20 CONTINUE 

IF(REF.GT.NPTS) GO TO 26 
STDDEV=SQRTF(SUMBDIF/(NPTS-1) ) 

GO TO 27 

26 STDDEV=SQRTF(SUMBDIF/NPTS) 

27 WRITE (61 *106) (E( I ) * 1=1*5) *STDDEV 
WRITE(61*111)RUN(REF) 

MNBDIFPC=SMBDIFPC/(NPTS-1 ) 

WRITE(61*112)MNBDIFPC 

GO TO 3 
51 REF=NPTS+1 

ALPHA ( REF ) =0 , 0 
VO ( REF ) =0 • 0 
LENGTH(REF)=0.0 
ROV I S ( REF ) =0 .0 
ROS I G ( REF ) =0 .0 
BFLASH(REF)=0.0 
DO 55 1=1, NPTS 

ALPHA(REF) = ALPHA (REF)+ALPHA( I ) 

VO(REF)=VO(REF)+VO( I > 

LENGTH ( REF )=LENGTH( REF )+LENGTH( I ) 

ROVIS(REF)=ROVlS(REF)+ROVIS( I ) 

ROS IG ( REF ) =ROSIG ( REF ) +ROS I G ( I ) 

BFLASH(REF)=BFLASH(REF)+BFLASH( I ) 

55 CONTINUE 

ALPHA ( REF )= ALPHA ( REF) /NPTS 
VO(REF)=VO(REF)/NPTS 
LENGTH ( REF) =LENGTH( REF) /NPTS 
ROV I S ( REF )=ROVIS( REF) /NPTS 
ROSIG ( REF )=ROSIG( REF ) /NPTS 
BFLASH ( REF )=BFLASH! REF) /NPTS 
RUN(REF)= 99999 
DPM ( REF ) = 1 • 

GO TO 24 

C IF IBRANCH EQUALS ZERO — DIMENSIONLESS NOS ARE RETURNED 

30 READ 107* IBRANCH 
I F ( E0F*60 )40,31 

31 IF ( IBRANCH. EQ. 0.0 ) GO TO 32 
GO TO 30 
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32 WR I TE ( 6 1 » 108 ) 

SUMKCMIN = 0».0 

THIS PROGRAM HAS THE OPTION OF USING REF RUN DATA CALCULATED IN 
THE BJORCK PROGRAM . 

CALC REF RUN DATA CARD (FROM BJORCK PROGRAM) SHOULD BE LAST DATA CARD IN 
DATA DECK 

IF t DPM ( REF ).EQ. 0.0 >34.33 

34 NPTS=NPTS-1 

33 DO 35 1=1 *NPTS 

KCM I N=2. 306658726*2. 0*32 • 172 5*DPM( I ) / ( DLlQ ( I ) *VO ( I )**2.0) 

SUMKCMI N=SUMKCM IN+KCMI N 
PRANDTL=1.0/(ALPHA( I )*ROVIS( I ) ) 

RED = 77.4192*ROVIS( I )*VO( I ) *DT S ( I ) 

REX =77 .4192*R0VIS( I ) *VO ( I ) *LENGTH ( I ) 

FROUDE=VO ( I )/SQRTF( LENGTH ( I ) *2. 68 104) 

WEBER=2359.737216*ROSIG( I ) *LENGTH ( I )*VO( I )**2.0 

peclet=prandtl*red 

WRITE (61 » 109) RUN ( I ) .ALPHA ( I ) » VO ( I ) » LENGTH ( I ) *ROVIS( I ) >ROSlG< I ) » 

1DTS ( I ) .KCMIN.PRANDTL.RED.REX.FROUDE, WEBER. PECLET 

35 CONTINUE 

avgkcmin=sumkcmin/npts 

WRITE (61 .110 ) AVGKCMIN 
GO TO 41 

10 1 FORMAT (A5»2F7.4,F8.3»E11.4,2F6.3*E11.4,F8,5*F8.4,2X»Il) 

102 FORM AT (5( F10.4) »6( 15) ) 

103 FORMAT!* El= * . F8 . 3 , 5X . *E2=* . F8 . 3 . 5X .*E3 = * » F8 . 3 . 5X .*E4=* . F8 . 3 . 

15X»*E5 = *»F8.3 .5X.*SUMBDIFSQ=*.E15.5 ) 

104 FORMAT ( *1RUN* ,5X.*ALPHA* »6X »*VO*.3X **LENGTH* » 6X »*ROV IS*»3X»*ROSIG* 

1 .3X»*BC0NV*»2X.*BFLASH*»4X»*BDIF*,2X.*BDIFPC*.5X.*SUMBDIFSQ*,6X. 
2*DTS*»9X. *NN* » / / ) 

105 FORMAT (A5.E11.4,F8.3,F7.3»E11.4,F8.5,F8.4,F8.4,F8.4»F8.1,E16.5, 

1 F8.4.9X.I1) 

106 FORMAT ( // / » 2X »*E1 = * » F8 . 3 » 5X . *E2=* » F8 . 3 » 5X »*E3=* . F8 . 3 . 5X »*E4=* . F8 . 3 
1.5X.*E5=* .F8.3.10X.*STDDEV=* »F10.4) 

107 FORMAT (II) 

10 8 FORMAT (*1RUN*»5X»*ALPHA*»6X.*V0*,3X.*LENGTH*.6X »*ROV IS*.3X,*ROSIG* 
1»5X.*DTS*»3X»*KCMIN*.4X»*PR*»6X»*RED*»10X.*REX*.8X.*FR*»9X.*WE*» 
27X»*PECLET*./ / ) 

109 FORMAT (A5»E11»4.F8»3.F7.3»E11»4.F8»5»F8»4.F8*3»F8.3>E12.2.E12.2»F 
18.2 .E12.2.E12.2 ) 

110 FORMAT </// .2X **AVGKCMI N = *»F8. 3.8 (/) ) 

111 FORMAT? ///,9X.*REF RUN = * » A5 » 2 ( / ) ) 

112 FORMAT ( 8X»*MNBDIFPC=*»F7.1 .8 ( / ) ) 

40 CALL EXIT 

END 

SCOPE 

• LOAD 

' RUN. 4, 10000.7 
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Bjorck Computer Listing 


' J0B.2750453.H0RD J.2 
< FTN»L.X»R 

PROGRAM BCONVECT 

DIMENSION RUN (405 ) »DPM(405) , DL IQ ( 405 ) » VO ( 405 ) ♦ ALPHA (405) *LENGTH 
1 (405) ,DTS(405) » R0VIS(405)» ROS I G ( 405 ) , BFLASH ( 405 ) »NN ( 405 ) , E ( 5 ) 

COMMON Ml»M»N»A(4l0*6)»B(4lC)»W(410),X(7) »RES(410) * V { £> * 6 ) »STD 
REAL KCMIN*LENGTH,MNBDIFPC 

INTEGER TERM1 , TERM2 ,TERM3 » TERM4 , TERM5 , WGHT , REF 
C THIS PROGRAM CORRELATES B-FLASH AND BCONVEC THEORIES— SIMILARITY 

C EQUATION— WITH EXPERIMENTAL DATA. 

C PROGRAM LOGI C9 (WITH BJORCK) DATA, EOF. NF CARD »NF CARD , NF CARD 

C ETC. , EOF, DATA, EOF, NF CARDS , EOF ,, EOF .( TAKES 2EOFS TOTERMINATE 

C PROGRAM.) FOR FORCE FIT ADD AN EXTRA CARD AFTER NF CARD — THIS CARD 

C WILL DEFINE THE TERM TO BE FIXED AND ITS EXPONENT VALUE. THIS CARD 

C IS REQUIRED ONLY IF Ml GREATER THAN O.FOR EXAMPLE TO FORCE-FIT THE 

C 1 »2ND»AND5TH TERMS WITH EXPONENTS OF .6,0.5. AND 0.4,— TYPE THE 

C TERM CARD AS FOLLOWS9 TERM 1=1, TERM 2=2 »TERM3=5 .LEAVE TERM4 AND 

C TERM5 BLANK AND TYPE E1=0 . 6 , E2=0 . 5 AND E3=0.4»LEAVE E4 ANDE5 BLANK 

C I • E • , THE E ( I ) MUST MATCH THE TERM BEING HELD CONSTANT. USE M1=0 

C FOR NO CONSTRAINTS AND THE TERM CARD IS NOT REQUIRED FOR A 3, 4, 

C OR 5 TERM FI T ( ALL EXPONENTS ARE JUGGLED). TO FORCE-FIT THE 3D TERM 

C WITH EXPONENTS. 5,SETM1 = 1 ON NF CARD THEN TYPE THE TERM CARD TERM1 

C =3 AND El=0 • 5 » ETC.. 

41 1 = 1 

1 READ 101 ,RUN( I ) ,DPM( I ) ,DLIQ( I ) ,VO( I ) ,ALPHA( I ) ,LENGTH( I ) ,DTS( I ) , 

IROVISt I ) , ROS I G ( I ) , BFLASH ( I ) , NN ( I ) 

IF (EOF, 60) 4, 2 

2 1 = 1+1 
GO TO 1 

4 NPTS= 1-1 

I F ( NPTS.EQ.O .0 ) GO TO 40 
M=NPTS 

THE FOLLOWING DO LOOP SHOULD BE ELIMINATED WHEN ROSIG IS TO BE USED 

THIS DO LOOP IS NECESSARY TO EVALUATE SCALE EFFECTS 

WEBER NO. WILL BE MEANINGLESS WHEN THIS DO LOOP IS USED. 

DO 44 I = 1 , NPTS 
ROS I G ( I ) = DTS ( I ) 

44 CONTINUE 

IF WGHT= 1 *W ( I ) = 1 AND IF WGHT=2» W ( I ) =BFLASH ( I ) *BFLASH ( I ) 

13 READ 110,NF,M1 , 1 BRANCH ,WGHT > REF 

IF REF =1 THE FIRST CARD READ WILL BE USED AS REF RUN, IF REF=10 THE 10TH 
CARD WILL BE USED, IF REF=0 THE REF RUN WILL BE CALCULATED AND REF WILL 
BE SET EQUAL TO NPTS+1. 

I F ( EOF, 60 ) 41 ,65 

65 IF(REF.EQ. 0.0)81, 5 

5 N=NF 

IF(Ml.GT.O.O) GO TO 45 

6 GO TO (66,76) , WGHT 

66 DO 8 1 = 1, M 

B ( I +M1 )=LOGF( BFLASH ( I ) /BFLASH ( REF ) > 

7 A ( I+Ml ,1 ) =LOGF( ALPHA (REF) /ALPHA ( I ) ) 

A ( I +M1 , 2 ) =LOGF ( VO ( I ) /VO (REF) ) 

A ( I+Ml ,3) =LOGF( LENGTH ( I ) /LENGTH (REF) ) 

A ( I +M1 ,4 ) =LOGF ( ROV I S ( I )/ROVIS(REF) ) 
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A ( I+Ml » 5 ) =LOGF ( ROSIG (I)/ROSIG(REF) ) 

8 CONTINUE 
CALL BJORCK 
PRINT 112 
12 DO 10 I*1»M 

JIM=I+M1 

10 PRINT 102»RUN(I) »RES(JIM) 

11 PRINT 103 . ( X ( J ) » J=1 *NF) 

PRINT 107 »STD 

14 PRINT 104 »NF 

C LOGF ( BCONVt I ) /BFLASH (REF) ) =ZZ 

SUMBDIF=0.0 
SMBDIFPC=0.0 
DO 18 1=1. M 

ZZ=0.0 

DO 16 J=1»NF 
16 ZZ=ZZ+A(I+M1.J)*X(J) 

BCONV=BFLASH ( REF) *EXPF ( ZZ ) 

BDIF=BCONV-BFLASH( I ) 

BDIFPC=(BDIF/BFLASH( I ) >*100.0 
SMBDIFPC=SMBDIFPC+ABSF(BDIFPC) 

SUMBDIF=SUMBDIF+BDIF*BDIF 

WRITE(61*105)RUN( I ) .ALPHA ( I ) .VO( I) . LENGTH ( I ) .ROVIS( I ) »ROSIG< I ) . 
1BCONV »BFLASH( I ) »BDIF » BD I FPC .SUMBD I F »DTS < I > .NN ( I > 

18 CONTINUE 

IF ( REF.GT .NPTS ) GO TO 20 
STDDEV=SQRTF(SUMBDIF/(NPTS-1 ) ) 

MNBDIFPC=SMBDIFPC/(NPTS-1) 

GO TO 21 

20 STDDEV=SQRTF ( SUMBD I F /NPTS ) 

MNBDIFPC=SMBDI FPC/ ( NPTS ) 

21 WRITEI61.106) STDDEV 

WRITE (61. 114) RUN (REF) .WGHT .MNBDI FPC 
C IF IBRANCH EQUALS ZERO — DIMENSIONLESS NOS ARE RETURNED 

I F ( IBRANCH.EQ.0.0 ) GO TO 32 
GO TO 13 

C A TERM AND EXPONENT CARD MUST FOLLOW ALL NF CARDS WHERE Ml .GT.O. 

45 READ 111.TERM1.TERM2.TERM3 .TERM4, TERM5 »(E( I) .1 = 1. Ml) 

I F ( EOF .60)40.47 
47 DO 60 K=1»M1 

GO TO (51.52.53.54.55) »K 

51 KK=TERM1 
GO TO 57 

52 KK=TERM2 
GO TO 57 

53 KK=TERM3 
GO TO 57 

54 KK=TERM4 
GO TO 57 

55 KK=TERM5 
GO TO 57 

57 B(K)=E(K) 

A(K.KK) =1.0 
KIM=KK-1 
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K I MN=KK+1 
DO 58 |_=1*KIM 

A ( K * L ) =0 « 0 

58 CONTINUE 
DO 59 L=KIMN*NF 
A ( K »L ) =0 • 0 

59 CONTINUE 

60 CONTINUE 
GO TO 6 

76 DO 78 1=1, M 

Bt I+Ml ) =LOGF( BFLASHt I ) /BFL ASH ( REF ) ). 

77 At I+Ml ,1 ) =LOGF( ALPHA (REF) /ALPHA ( I ) ) 

At I+Ml ,2) =LOGF(VO( I) /VO (REF) ) 

At I+Ml, 3 )=LOGF( LENGTH ( I ) /LENGTH ( REF ) ) 

At I+Ml >4) =LOGF( ROVIS ( I ) /ROVIS(REF) > 

At I+Ml *5) =LOGF(ROSIG{ I ) /ROSIG(REF) ) 

W ( I+Ml ) =THE WEIGHT REQUIRED FOR LEAST SQUARES MINIMIZATION ON THE 

LINEAR PROBLEM NOT MINIMIZATION IN THE LOGARITHMIC SENSE 

MINIMIZING IN THE LOG SENSE IS DONE ON BD I FPC*BD I FPC ( WHEN WGHT=1) 

W( I+Ml ) =BFLASH ( I ) +BFLASH ( I ) 

78 CONTINUE 
CALL WEIGHT 
PRINT 112 
GO TO 12 

C STATEMENT 81 STARTS CALC OF REF CONDITIONS USING BASIC BCONV EQUATION 

81 ISETM1=M1 
Ml = 0 « 

N=6 

DO 82 I = 1 , NPTS 

Bt I ) =LOGF ( BFLASH ( I ) ) 

At I ,1 ) = 1.0 

At I * 2 ) =LOGF ( ALPHA ( I ) ) 

At I * 3 ) =LOGF ( VO ( I ) ) 

At I ,4)=L0GF(LENGTH( I ) ) 

At I ,5 )=LOGF( ROVISt I ) ) 

At I , 6 ) =LOGF ( ROS I G ( I ) ) 

82 CONTINUE 
CALL BJORCK 

83 REF=NPTS+1 

ALPHA ( REF ) =0 • 0 
VO(REF) =0.0 
LENGTH ( REF ) =0 . 0 
ROV I S ( REF ) =0 . 0 
ROS I G ( REF ) =0 . 0 
DO 85 1=1, NPTS 

ALPHA(REF)=ALPHA(REF)+ALPHA( I ) 

VO(REF)=VO(REF)+VO( I ) 

LENGTH ( REF )=LENGTH( REF )+LENGTH( I ) 

ROVISt REF) =ROVI St REF l+ROVI St I ) 

ROSIG(REF)=ROSIG( REF)+ROSIG( I) 

85 CONTINUE 

ALPHA ( REF )=ALPHA( REF) /NPTS 
VOt REF) =V0( REF) /NPTS 
LENGTH! REF) =LENGTH( REF) /NPTS 
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ROVI S ( REF )=ROVIS( REF) /NPTS 
ROS IG ( REF ) =ROS IG(REF)/NPTS 

PP = X ( 1 ) +X ( 2 ) *LOGF ( ALPHA ( REF ) ) +X ( 3 ) *LOGF ( VO ( REF ) ) +X ( 4 ) *LOGF ( LENGTH 
1 (REF) )+X(5>*LOGF(ROVIS(REF) ) +X ( 6 > *LOGF ( ROS I G ( REF ) ) 
BFLASH(REF)=EXPF(PP) 

M1=ISETM1 
RUN (REF)=99999 
DPM ( REF ) =0. 0 
DLIQ( REF) =0.0 
DTS ( REF ) =DTS ( 1 ) 

NN(REF)=NN( 1) 

PRINT 115 »BFLASH( REF) » ALPHA ( REF) * VO (REF) * LENGTH (REF) .ROVIS(REF) * 

1 ROSIG(REF) 

A PUNCH STATEMENT IDENTICAL TO PRINT 115 SUPPLIES THIS CALCULATED 
REF DATA TO VAR I AT I ON-OF-PA-RAMETER PROGRAM 
GO TO 5 

32 WRITE(61»108) 

SUMKCMIN=0.0 

33 DO 35 I = 1 ♦ NPTS 

KCMIN=2.306658726*2.0*32.1725*DPM( I )/(DLlQ( I )*VO( I )**2.0) 

SUMKCM IN=SUMKCM IN+KCMI N 
PRANDTL=1.0/ (ALPHA ( I ) *ROVIS ( I ) ) 

RED=77.4192*ROVIS( I ) *VO ( I )*DTS( I ) 

REX = 77 .4192*ROVI S ( I ) *VO ( I ) *LENGTH ( I ) 

FROUDE=VO ( I ) /SQRTF( LENGTH ( I ) *2.68104) 

WEBER=2359.737216*ROSIG( I ) ^LENGTH ( I ) *VO ( I ) **2.0 
PECLET=PRANDTL*RED 

WRITE ( 61 .109) RUN ( I ) , ALPHA ( I ) » VO ( I ) .LENGTH ( I ) ,ROVIS( I ) »ROSlG( I ) . 
1DTS ( I ) > KCMI N » PRANDTL.RED * REX >FROUDE> WEBER. PECLET 
35 CONTINUE 

AVGKCM I N= SUMKCM IN /NPTS 
WRITE (61 .113) AVGKCM IN 

101 FORMAT! A5 .2F7.4.F8.3 » E 1 1 . 4 , 2F6 . 3 , El 1 .4 . F8 . 5 * F8 . 4 . 2X . 1 1 ) 

102 FORMAT (A5 .E20.5 ) 

103 FORMAT(*l ZZ IS THE SUM OF THE FOLLOWING TERMS9*//*0 ZZ=*El7.10. 

1 8H*A ( I » 1 ) E 17. 10 .8H*A ( I . 2 ) E17. 10 »8H*A ( I . 3) 

2 /*0*E17»10.8H*A( I. 4 ) E 1 7 . 10 , 8H*A ( I » 5 ) E 1 7 . 10 , 8H* A < I * 6)) 

104 FORMAT(*l*» II ,* FUNCTIONS GIVES THE FOLLOWING FIT9*///*RUN*»5X. 
l*ALPHA*»6X»*VO*»3X.*LENGTH*.6X.*ROVI5*»3X»*ROSIG*»3X»*BCONV*,2X> 

2 *BFLASH*.4X,*BDIF*»2X.*BDIFPC*>5X»*5UMBDIFSQ*»6X,*DTS*.9X,*NN*//) 

105 FORMAT (A5»E11.4*F8»3»F7»3»E11»4*F8»5»F8.4*F8.4*F8.4»F8«1»E16.5» 

1 F8.4.9X* 1 1 > 

106 FORMAT ( ////.*STD DEVIATION AS COMPUTED BY SUMBDI FSQ=* »F10. 5 ) 

107 FORMAT ( *0STD DEV OF THE FIT=*.F10.5) 

108 FORMAT ( *1RUN*,5X» *ALPHA* »6X » *VO* »3X»*LENGTH*.6X »*ROV IS*»3X»*ROSIG* 

♦6X,*RED*.lOX»*REX*»8X»*FR*.9X.*WE*. 

1 .4.F8.5 ,F8.4,F8.3,F8.3»E12.2»E12.2»F 


.3) 

» 20X ♦ *WGHT =* » 15 .20X»*MNBDIFPC = *»F7.1 
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1>5X.*DTS*»3X»*KCMIN*.4X.*PR* 
27X.*PECLET*.// ) 

109 F0RMAT(A5.E11»4»F8.3»F7»3.E1 
18.2 .E12.2 .E12.2 ) 

110 FORMAT ( 5 ( I5> ) 

111 FORMAT (5(I5)»5(F10.3)) 

112 F0RMAT(*1RUN*,15X,*RES*.//) 

113 FORMAT ( ///,2X »*AVGKCMIN=* >F8 

114 FORMAT (///.2X»* REF RUN=*.A5 
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1 f 6 ( / ) ) 

115 FORMAT ( *CALC BFLASH ( REF ) =* > F8 .4 * / f *CALC ALPHA ( REF ) =* , El 1 .4 , / , *CALC 
1V0(REF)=*,F8.3>/**CALC LENGTH ( REF J = *, F7. 3 ,/ ,*CALC ROV I S ( REF ) =* , 


2 El 1 • 4 » / * *CALC ROS I G ( REF ) =* , F8 • 5 , 6 ( / ) ) 

40 CALL EXIT 
END 

SUBROUTINE BJORCK BJOROOl 

DIMENSION Q(410,6),G(6),F(410),Y(6),D(6) ,NPV(6) , P ( 7 ) , SW ( 4 10 ) 

COMMON Ml ,M » N , A ( 410 » 6 ) >B(410),W(410),X(7) ,RES(4lO) ,V(6,6) >STD 


BJOR004 
BJOR005 
BJOR006 
BJOR007 
BJ0R008 
BJOR009 
BJOROIO 
BJOROll 
BJOR012 
BJOR013 
BJOR014 
BJOR015 
BJOR016 
BJ0R017 
BJOR018 
BJOR019 
BJOR020 
BJOR021 
BJOR022 
BJ0R023 
BJOR024 
BJOR025 
BJOR026 

THE SUBROUTINE ACPROD ( A , B ,N » C ) SUMS THE FIRST N PRODUCTS OF BJOR027 
A ( I ) *B ( I ) AND ADDS C IN DOUBLE THE OPERATING PRECISION AND IS»BJOR028 


THEREFORE, MACHINE DEPENDENT. THE OPERATING PRECISION IS BJOR029 

DEFINED BY THE VALUE OF ETA WHICH IS 1/2 RAISED TO A POWER BJOR030 
EQUAL TO THE NUMBER OF BITS CARRIED IN THE ARITHMETIC BJOR031 

OPERATIONS (IN THIS CASE, 36) BJOR032 

BJOR033 

DATA (ETA=1.45519152284E“11 ) BJOR034 

DO 10 1=1, M BJOR035 

10 SW( I ) =1. BJOR036 

11 ETA2=ETA*ETA BJOR037 


THE FOLLOWING VALUES MUST BE SET BEFORE THE CALL 

A ( I , J ) = THE J FUNCTION EVALUATED AT THE I POINT 

B ( I ) = THE VALUE OF THE I ORDINATE 

W(I> = THE WEIGHT TO BE APPLIED TO THE I POINT 

M = THE NUMBER OF POINTS ( I ) 

N = THE NUMBER OF FUNCTIONS ( J ) 

Ml = THE NUMBER OF ROWS IN A THAT MUST BE SATISFIED EXACTLY 
(NUMBER OF CONSTRAINTS) 

AFTER THE CALL 

X(J) = THE MEAN- LEAST-SQUARE COEFFICIENT OF THE J FUNCTION 
RES ( I ) = THE MEAN LEAST SQUARE DEVIATION OF THE I POINT 
V(I,J) = THE VARIANCE - COVARIANCE MATRIX 
STD = THE STANDARD DEVIATION OF THE FIT 

THE METHOD IS THAT OF AKE BJORCK, AS DESCRIBED IN 
NORDISK TIDSKRIFT FOR I NFORMAT I ONS-BEHANDL I NG 
BIT VOL. 7, PP. 257-278 (1967) AND V0L#8, PP.8-30 (1968) 

THE WEIGHTING OF EACH POINT IS 1.0 UNLESS THE ENTRY WEIGHT 
IS USED. 


C = 36*N BJOR038 

TOL2=C*ETA2 BJOR039 

MP=M1+1 BJOR040 

DO 13 J= 1 »N BJOR041 

NPV ( J ) = J BJOR042 

DO 12 1 = 1, N BJOR043 

12 V( I , J ) =0 • BJOR044 

DO 13 1=1, M BJOR045 

13 Q( I 9 J ) =A ( I , J ) *SW ( I ) BJOR046 

MV= 1 BJOR047 

MH=M 1 BJOR048 
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NTR= 1 

DO 30 K= 1 * N 
NP = K 

I F ( K • NE*MP ) GO TO 15 
M V = MP 
MH = M 

14 NTR = 1 

15 DS=0« 

DO 20 J=K>N 
L = NPV ( J ) 

I F ( NTR «EQ#0 ) GO TO 17 
D ( J ) =0 • 

DO 16 I =MV*MH 

16 D( J )=D( J) +Q( I »L)*Q( I »L) 

17 I F ( DS*GE«D ( J ) ) GO TO 20 
DS=D ( J ) 

NP = J 

20 CONTINUE 
L=NPV ( NP ) 

IF(NTR.EQ.l) dm=ds 
NTR = 0 

I F { DS • LT • ETA*DM ) GO TO 14 
IF(NP.EQtK) GO TO 22 
D ( NP ) =D ( K ) 

NPV ( NP > = NPV ( K ) 

NPV ( K ) =L 
KM=K-1 

DO 21 1=1 »KM 
C=V ( I *NP) 

V( I »NP)=V( I >K) 

21 V ( I *K ) =C 

22 DS=0. 

DO 23 I=MV,MH 

23 DS=DS+Q ( I *L ) *Q ( I*L) 

D ( K ) =DS 

IF(K.LE.Ml) GO TO 24 

IF< DS.LE»T0L2*D(MP)*V(MP»K)*V(MP>K) ) GO TO 31 

24 IF(DS.LE.O.) GO TO 44 
KM=K+1 

IF(KM.GT.N) GO TO 30 
DO 27 J = KM * N 
NP=NPV( J) 

DO 25 I=MV,MH 

25 V(K*J)=V(K*J)+Q( I»L)*Q( IfNP) 

V<K*J)=V<K»J)/DS 

DO 26 I=MV*M 

26 Q( I *NP)=Q( I tNP)-V(Kf J)*Q< I »L) 

27 D ( J ) = D ( J ) *-DS*V 

30 N1=K 

31 DS= 1 • 

33 STD=0. 

DO 35 I = 1 *M 
RES < I ) = 0 • 

35 F ( I ) =B ( I ) *SW { I ) 


BJ0R049 

BJOR050 

BJOR051 

BJOR052 

BJOR053 

BJ0R054 

BJOR055 

BJOR056 

BJOR057 

BJ0R058 

BJOR059 

BJ0R060 

BJOR061 

BJ0R062 

BJOR063 

BJOR064 

BJOR065 

BJOR066 

BJOR067 

BJOR068 

BJOR069 

BJOR070 

BJOR071 

BJ0R072 

BJOR073 

BJ0R074 

BJOR075 

BJOR076 

BJOR077 

BJOR078 

BJOR079 

BJOR080 

BJ0R081 

BJOR082 

BJOR083 

BJOR084 

BJOR085 

BJOR086 

BJOR087 

BJOR088 

BJ0R089 

BJOR090 

BJOR091 

BJOR092 

BJOR093 

BJOR094 

BJOR095 

BJOR096 

BJ0R097 

BJOR098 

BJOR099 

BJOR100 

BJORIOI 

BJOR102 
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79 


L 




DO 37 J=1 »N 

BJOR103 


X( J>=0. 

BJOR104 

37 

G(J)=0. 

BJOR105 


X ( N+l ) =-l • 

BJOR106 


I TER=-1 

BJOR107 


IF(DS.GT.O.) GO TO 55 

BJOR108 


DO 39 J=1*N 

BJOR109 


DO 39 1=1 »M1 

BJORllO 

39 

V( I »J)=0. 

BJ0R111 


RETURN 

BJ0R112 

40 

FORMAT ( 12OH0$$$$$$ AT LEAST TWO OF THE CONSTRAINTS ARE EQUI VALENT. BJOR113 


1 FITTING WAS TERMINATED 

AND ALL COEFFICIENTS SET TO ZERO $$$$$$$) BJ0R1 14 

41 

FORMAT ( 120H0 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ ITERATION FAILBJORH5 


IS TO IMPROVE THE SOLUTION $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$) BJOR116 


ENTRY WEIGHT 

BJ0R117 


DO 43 I = 1 »M 

BJ0R118 


I F ( W { I ) .LT.O. ) W(I)=0. 

BJ0R119 

43 

SW( I ) =SQRT ( W ( I ) ) 

BJOR120 


GO TO 11 

BJ0R121 

44 

PRINT 40 

BJ0R122 


GO TO 33 

BJ0R123 

45 

DO 46 I = 1 »M 

BJ0R124 

46 

RES ( I ) =RES ( I )+F( I) 

BJOR125 


DO 47 K=1 »N1 

BJ0R126 


J=NPV ( K ) 

BJ0R127 

47 

X(J)=X< J)+G(K) 

BJOR128 


IF(DX2.LT.ETA2*X1 ) GO TO 

70 BJ0R129 


IF ( DR2.LT »ETA2*R1) GO TO 

70 BJOR130 


IF(64.D0*DX2.LT.DX1) GO 

TO 50 BJ0R131 


IF(64.D0*DR2.LT.DR1 > GO 

TO 50 BJOR132 


PRINT 41 

BJOR133 


GO TO 70 

BJ0R134 

50 

DX1=DX2 

BJOR135 


DR1=DR2 

BJ0R136 


DO 52 K=1 »N1 

BJ0R137 


J=NPV ( K ) 

BJOR138 


DO 51 I = 1 »M 

BJ0R139 

51 

F ( I ) =A ( I ♦ J)*SW( I > 

BJOR140 

52 

G ( K ) =-ACPROD ( F » RES »M »0 • ) 

BJ0R141 


C = 0. 

BJ0R142 


DO 54 I = 1 *M 

BJ0R143 


DO 53 J=1 »N 

BJ0R144 

53 

P(J)=A( I*J)*SW( I) 

BJ0R145 


P ( N+l ) =B( I )*SW( I ) 

BJ0R146 


IF(I.GT.Ml) C = RES ( I ) 

BJ0R147 

54 

F ( I ) = -AC PROD (P»X»N+1 * C ) 

BJ0R148 

55 

MV=1 

BJ0R149 


MH=M1 

BJOR150 


DO 60 K=l»Nl 

BJ0R151 


J=NPV ( K ) 

BJOR152 


IF(K.NE.MP) GO TO 56 

BJOR153 


MV=MP 

BJ0R154 


MH = M 

BJOR155 

56 

Y { K ) =G ( K ) 

BJOR156 

08 

05 



80 




IF(K.EQ.l) GO TO 58 

BJOR157 


L = K-1 

BJOR158 


DO 57 1=1 tL 

BJ0R159 

57 

Y<K)=Y(K)-V( I.K)*Y( I ) 

BJOR160 

58 

G ( K ) =0 • 

BJ0R161 


IF(K.GT.Ml) G(K)=-Y(K) 

BJOR162 


DO 59 I =MV *MH 

BJOR163 

59 

G ( K ) =G t K ) +0 ( I » J ) *F f I ) 

8J0R164 


G(K)=G<K)/D<K) 

BJ0R165 


DO 60 I =MV *M 

BJ0R166 

60 

F < I ) =F ( I )-G( K ) *Q( I ,J) 

BJ0R167 


IF(Ml.LE.O) GO TO 65 

BJ0R168 


DO 62 1=1 »M1 

BJ0R169 

62 

F ( I ) =0* 

BJOR170 


DO 64 K=1»M1 

BJ0R171 


J=NPV(<) 

BJ0R172 


C=0. 

BJ0R173 


DO 63 I = 1 *M 

BJ0R174 

63 

C=C+Q ( I * J)*F( I ) 

BJOR175 


C= ( C-Y ( K ) )/D(K) 

BJOR176 


DO 64 1=1 *M1 

BJ0R177 

64 

F ( I ) =F ( I )-C*Q( I »J) 

BJOR178 

65 

DO 67 LL = 2 * N1 

BJOR179 


K=N1-LL+1 

BJOR180 


L = K+1 

BJ0R181 


DO 66 I=L»N1 

BJOR182 

66 

G ( K ) =G ( K ) -V ( K > I )*G( I ) 

BJOR183 

67 

CONTINUE 

BJOR184 


DX2=0. 

BJOR185 


DO 68 K=1 »Nl 

BJ0R186 

6 8 

DX2=DX2+G(K)*G(K) 

BJOR187 


DR2 =0 • 

BJOR188 


DO 69 I = 1 *M 

BJ0R189 

69 

DR2=DR2+F ( I ) *F ( I ) 

BJOR190 


ITER=ITER+1 

BJ0R191 


I F t ITER.GT.O) GO TO 45 

BJOR192 


X1=DX2 

BJOR193 


R1=DR2 

BJ0R194 


DX 1 =X 1*65 • 

BJ0R195 


GO TO 45 

BJ0R196 

70 

NRP = 0 

BJ0R197 


DO 73 I - 1 »M 

BJOR198 


IF(I.LE.Ml) SW C I ) =0 . 

BJOR199 


IFfSWm « GT • 0 • ) GO TO 72 

BJOR200 


DO 71 J= 1 *N 

BJOR201 

71 

P( J>=A( I.J) 

BJOR202 


P(N+1)=B( I ) 

BJOR203 


RES ( I ) =-ACPROD( P*X»N+1 ,0. ) 

BJOR204 


GO TO 73 

BJOR205 

72 

NRP=NRP+1 

BJOR206 


STD=STD+RES< I )*RES( I ) 

BJ0R207 


RES ( I >=RES( I ) /SW( I ) 

BJ0R208 

73 

CONTINUE 

BJOR209 


I =M1 + NRP-N 1 

BJOR210 


71 08 05 


81 



non 


IF(I.GT.O) GO TO 74 BJ0R211 

DS=0. BJ0R212 

GO TO 75 BJOR213 

74 DS=STD/I BJOR214 

75 STD=SQRT(DS) BJOR215 

L=N1 BJOR216 

IF(L*GT.M1) V(L*L)=1./D(L) BJOR217 

I F ( N1 • GE • N ) GO TO 77 BJOR218 

DO 76 I=N 1 *N BJOR219 

76 V ( L * I ) “0 • BJOR220 


77 


78 


79 

80 

81 


82 

83 


84 

85 


89 


J = L-1 

IF(J.GT.Ml) V(J,J)=1./D< J) 

DO 78 K = L *N 1 
P(K)=V( J,K) 

I =N 1 

DO 80 KA = J * N 
DO 79 K=L *N 

V(ItJ)=V( ItJ)-P(K)*V(KtI) 

1 = 1-1 

DO 81 K=L »N 
V( J*K)=V(K>J) 

L=L-1 

IF(L.GT.l) GO TO 77 
DO 83 I = 1 »N 
DO 82 J= 1 *N 
K=NPV( J) 

P(K)=V( It J) 

DO 83 J= 1 *N 
V ( I » J ) =P ( J ) 

DO 85 I = 1 »N 
DO 84 J= 1 *N 
K=NPV ( J ) 

P(K)=V( Jtl) 

DO 85 J= 1 *N 
V ( J 1 1 ) =P( J ) 

DO 89 J= 1 »N 

DO 89 1=1 tN 

V( I ♦J)=V( I» J)*DS 

RETURN 

END 

FUNCTION AC PROD ( AA >BB *N * CC ) 

THIS SUBROUTINE FORMS THE DOUBLE 
FIRST N VALUES OF THE ARRAYS » AA 


BJOR221 

BJOR222 

BJOR223 

BJOR224 

BJOR225 

BJOR226 

BJOR227 

BJOR228 

BJOR229 

BJOR230 

BJOR231 

BJOR232 

BJOR233 

BJOR234 

BJOR235 

BJOR236 

BJOR237 

BJOR238 

BJOR239 

BJOR240 

BJOR241 

BJOR242 

BJOR243 

BJOR244 

BJOR245 

BJOR246 

BJOR247 

BJOR248 

BJOR249 

BJOR250 

BJOR251 

PRECISION PRODUCT OF THE BJOR252 

AND BB. CC IS THEN ADDED IN BJOR253 


DOUBLE PRECISION AND THE RESULT 
DIMENSION AA ( 9 ) *BB(9) 

DOUBLE PRECISION StAtB*C 
S=0 #D0 
C = CC 

DO 1 I = 1 ♦ N 
A = AA ( I ) 

B = BB ( I ) 

1 S=S+A*B 
ACPROD=S+C 
RETURN 


RETURNED IN SINGLE PRECISION BJOR254 

BJOR255 

BJOR256 

BJOR257 

BJOR258 

BJOR259 

BJOR260 

BJOR261 

BJOR262 

BJOR263 

BJOR264 
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END 


BJOR265 


SCOPE 

' LOAD 

• RUN » 2 ♦ 10000 * 7 
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